#include <fcntl.h>
#include <ctype.h>
#include  "mcio.h" 

int  usecm = NO ;
int plinkinputmode = NO ;
static int snprawtab = NO ;
static int debug = NO ;
extern  char *trashdir ;
extern int qtmode ;
extern int verbose ;
extern int familynames ;
extern double lp1, lp2  ;
extern double a1, b1  ;

extern int packmode ;  
extern char *packgenos ;
extern char *packepath ;
extern long packlen, rlen ;
extern int malexhet ;
extern int hashcheck ;
extern int outputall ;
extern int sevencolumnped ;
static int dofreeped = YES ;

int tempnum = 0 ;
int tempfake = 0 ;

static int *pedcols = NULL ;
static int numpedcols = 0 ;

static int *pedcolsa[3] ;  
static int numpedcolsa[3] ;  


static int maxgenolinelength = -1 ;
static int tersemode = NO ;
int checksizemode = YES ;
int pedignore = YES ;
enum outputmodetype outputmode = PACKEDANCESTRYMAP ;
static double maxgpos[MAXCH] ;

int getbedgenos(char *gname, SNP **snpmarkers, Indiv **indivmarkers, 
 int numsnps, int numindivs, int nignore)  ;

void freeped()  ;

static char x2base(int x) ;
static void gtox(int g, char *cvals, int *p1, int *p2)  ;

int ancval(int x)  ;
static int setskipit(char *sx) ;  // ignore lines in snp, map files
int calcishash(SNP **snpm, Indiv **indiv, int numsnps, int numind, int *pihash, int *pshash) ;  

SNPDATA  *tsdpt ;

int getsnps(char *snpfname, SNP ***snpmarkpt, double spacing,
  char *badsnpname, int *numignore, int numrisks) 
// returns number of SNPS
// numrisks 
{
/* read file of real SNPS store in temporary structure */

 SNPDATA **snpraw, *sdpt ;
 static SNP **snpmarkers ;
 SNP *cupt ;
 int **snppos ;
 int nreal, nfake, numsnps = 0, i, t ;
 int *snpindx ;
 double xspace ;
 char sss[20] ;

 if (snpfname == NULL) fatalx("(getsnps) null snpname") ;
 xspace = spacing ;
 nreal = getsizex(snpfname) ;
 if (nreal <= 0) fatalx("no snps found: snpfname: %s\n", snpfname) ;
 ZALLOC(snpraw, nreal, SNPDATA *) ;
 if (pedcols == NULL) {
  ZALLOC(pedcols, nreal, int) ;  
  ivclear(pedcols, -1, nreal) ;
  numpedcols = nreal ;
 }
 for (i=0; i<nreal ; i++) { 
  ZALLOC(snpraw[i], 1, SNPDATA) ;
  snpraw[i] -> inputrow = -1 ;
  snpraw[i] -> alleles[0] = '1' ;
  snpraw[i] -> alleles[1] = '2' ;
 }
 readsnpdata(snpraw, snpfname) ;
 dobadsnps(snpraw, nreal, badsnpname) ;

 ZALLOC(snppos, nreal, int *) ;
 for (i=0; i<nreal; i++) {  
  ZALLOC(snppos[i], 3, int) ;
 }


 for (i=0; i<nreal ; i++) {
  sdpt = snpraw[i] ;
  snppos[i][0] = sdpt -> chrom ;   
//  if (sdpt->ignore) snppos[i][0] = 99 ;
  sprintf(sss, "%15.6f", sdpt -> gpos) ;
  sscanf(sss, "15.6f", &sdpt -> gpos)  ;
  t = snppos[i][1] = nnint((sdpt -> gpos)*GDISMUL) ;       
  snppos[i][2] = nnint(sdpt -> ppos) ;      
// sdpt -> gpos = ((double) t)/ GDISMUL ;
 }

 ZALLOC(snpindx, nreal, int) ;
 ipsortit(snppos, snpindx, nreal, 3) ;

 if ((usecm) && (xspace>0.5)) { 
   printf("*** warning fake spacing given in cM\n") ;
   xspace /= 100.0 ;
 }
 
// get number of fakes

 nfake = numfakes(snpraw, snpindx, nreal, xspace) ;
 numsnps = nreal + nfake ;
 
 tempnum = numsnps ;
 tempfake = nfake ;

// allocate storage

 ZALLOC(snpmarkers, numsnps, SNP *) ;
 for (i=0; i<numsnps; i++) {  
  ZALLOC(snpmarkers[i], 1, SNP) ;
  cupt = snpmarkers[i] ;
  clearsnp(cupt) ;
  ZALLOC(cupt -> modelscores, numrisks, double) ;
  ZALLOC(cupt -> totmodelscores, numrisks, double) ;
 }
 tsdpt = snpraw[0] ;
 *snpmarkpt = snpmarkers ;
  numsnps = loadsnps(snpmarkers, snpraw, snpindx, nreal, xspace, numignore) ;


// and free up temporary storage
 for (i=0; i<nreal ; i++) { 
  free(snpraw[i]) ;             
  free(snppos[i]) ;             
 }
 free(snpraw) ;
 free(snppos) ;
 free(snpindx);

 /* printf("numsnps: %d\n", numsnps) ; */

/**
 if (pedcols != NULL) { 
  printimat(pedcols, 1, MIN(100, numsnps)) ;
 }
*/
 cupt = snpmarkers[0] ; 
 if (isnumword(cupt -> ID)) printf("*** warning: first snp is number.  perhaps you are using .map format\n") ;

 return numsnps ;
}


int getsizex(char *fname) 
// number of real lines 
{
 char line[MAXSTR+1], c ;
 char *spt[MAXFF], *sx ;
 int nsplit, num=0 ;
 int skipit ;
 int len ;

 FILE *fff ;
 openit(fname, &fff, "r") ;
 line[MAXSTR] = '\0' ;
 while (fgets(line, MAXSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit == 0) continue ;
   sx = spt[0] ;
   skipit =  setskipit(sx) ;
   if (skipit == NO) {
    ++num ;
   }
// now flush the rest of the line if necessary.
   len = strlen(line) ;  
   c = line[len-1] ;
   if (c != '\n') { 
     while ((c = fgetc(fff)) != EOF) {  
      if (c == '\n') break ;
     }
   }
   freeup(spt, nsplit) ;
   continue ;
 }
 fclose(fff) ;
 fflush(stdout) ;
 return num ;
}

int ismapfile(char *fname) 
// PLINK map file ? 
// just look at file name (perhaps should look at format) 
{

 char *sx ;
 int len ;
 len = strlen(fname) ;
 if (len<4) return NO ;
 sx = fname+len-4 ;

 if (strcmp(sx, ".map") == 0) return YES  ;
 if (strcmp(sx, ".bim") == 0) return YES  ;
 
 if (len<7) return NO ;
 sx = fname+len-7 ;
 if (strcmp(sx, ".pedsnp") == 0) return YES  ;

 return NO ;

}
int ispedfile(char *fname) 
// PLINK ped file ? 
// just look at file name (perhaps should look at format) 
{

 char *sx ;
 int len ;
 len = strlen(fname) ;
 if (len<4) return NO ;
 sx = fname+len-4 ;

 if (strcmp(sx, ".ped") == 0) return YES  ;
 if (strcmp(sx, ".fam") == 0) return YES  ;

 if (len<7) return NO ;
 sx = fname+len-7 ;
 if (strcmp(sx, ".pedind") == 0) return YES  ;

 return NO ;

}

int isbedfile(char *fname) 
// PLINK ped file ? 
// just look at file name (perhaps should look at format) 
{

 char *sx ;
 int len ;
 len = strlen(fname) ;
 if (len<4) return NO ;
 sx = fname+len-4 ;

 if (strcmp(sx, ".bed") == 0) return YES  ;
 return NO ;

}

int readsnpdata(SNPDATA **snpraw, char *fname) 
{
 char line[LONGSTR] ;
 char *spt[MAXFF], *sx ;
 int nsplit, num=0, k ;
 int skipit ;
 SNPDATA *sdpt ;

 double maxg = -9999.0 ;

 FILE *fff ;
 int chrom ;
 int nbad = 0 ;

 if (ismapfile (fname)) {  
    plinkinputmode = YES ;
    return readsnpmapdata(snpraw, fname)  ;
 }
 

 vclear(maxgpos, -9999.0, MAXCH) ; 
 openit(fname, &fff, "r") ;
 while (fgets(line, LONGSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit == 0) continue ;
   sx = spt[0] ;
   skipit = setskipit(sx) ;
   if (skipit == NO) { 
    if (nsplit<4) fatalx("(readsnpdata) bad line: %s\n",line) ;
    sdpt = snpraw[num] ;
    sdpt -> inputrow = num ;

    if (strlen(spt[0]) >= IDSIZE) fatalx("ID too long\n", spt[0]) ;
    strcpy(sdpt->ID, spt[0]) ;

    
    sdpt -> chrom = chrom = str2chrom(spt[1]) ; 

    if ((chrom>=MAXCH)  || (chrom <=0)) {
     if (nbad<10) printf("warning: bad chrom: %s", line) ;
     ++nbad ; 

     sdpt -> chrom = MIN(chrom, BADCHROM) ;  
     sdpt -> chrom = MAX(chrom, 0) ;  
     
     sdpt -> ignore = YES ;  
    }

    sdpt -> gpos = atof(spt[2]) ;
    if (sdpt->gpos > 100) {
      if (sdpt->gpos > 1.0e6) 
        fatalx("absurd genetic distance:\n%s\n", line) ;
     if (!usecm) {
      printf("*** warning.  genetic distances are in cM not Morgans\n") ;
      printf("%s\n",line) ;
     }
     usecm = YES ;
    }
    maxgpos[chrom] = MAX(maxgpos[chrom], sdpt -> gpos) ;
    maxg = MAX(maxg, maxgpos[chrom]) ;

    sdpt -> ppos = atof(spt[3]) ;
    if (nsplit<8) { 
     ivzero(sdpt->nn,4) ; 
     if (nsplit==6) {  
      sx = spt[4] ; sdpt -> alleles[0] = sx[0] ;
      sx = spt[5] ; sdpt -> alleles[1] = sx[0] ;
     }
    }
    else {
     for (k=0; k<4; k++) { 
      sdpt->nn[k] = atoi(spt[4+k]) ;
     }
     if (nsplit==10) {  
      sx = spt[8] ; sdpt -> alleles[0] = sx[0] ;
      sx = spt[9] ; sdpt -> alleles[1] = sx[0] ;
     }
    }
    ++num ;
   }
   freeup(spt, nsplit) ;
   continue ;
 }

 if (maxg<=0.00001)  {
  printf("genetic distance set from physical distance\n") ;
  usecm = NO ;
  for (k=0; k<num ; ++k) {
   snpraw[k] -> gpos = 1.0e-8 * snpraw[k] -> ppos ;
  }
 }

 if (usecm) {
  for (k=0; k<num ; ++k) {
   snpraw[k] -> gpos /= 100.0 ;
  }
 }
 
 fclose(fff) ;
 return num ;
}

int readsnpmapdata(SNPDATA **snpraw, char *fname) 
{
 char line[MAXSTR] ;
 char *spt[MAXFF], *sx ;
 int nsplit, num=0, k, t ;
 int skipit, len ;
 SNPDATA *sdpt ;
 int nbad = 0 ;

 FILE *fff ;
 int chrom ;
 double maxg = -9999.0 ;


 vclear(maxgpos, -9999.0, MAXCH) ; 
 openit(fname, &fff, "r") ;
 while (fgets(line, MAXSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit == 0) continue ;
   sx = spt[0] ;
   skipit = setskipit(sx) ;
   if (skipit == NO) { 
    if (nsplit<4) fatalx("(readsnpmapdata) bad line: %s\n",line) ;
    sdpt = snpraw[num] ;
    if (strlen(spt[1]) >= IDSIZE) fatalx("ID too long\n", spt[1]) ;
    strcpy(sdpt->ID, spt[1]) ;
   

    if (nsplit >=6 ) {  
     sx = spt[4] ; sdpt -> alleles[0] = sx[0] ;
     sx = spt[5] ; sdpt -> alleles[1] = sx[0] ;
     if (sdpt->alleles[0] == '0') sdpt -> alleles[0] = 'X' ; 
     if (sdpt->alleles[1] == '0') sdpt -> alleles[1] = 'X' ; 
    }
    else { 
     cclear((unsigned char *) sdpt -> alleles, CNULL,  2) ;
    }


    sx = spt[0] ;
    sdpt -> chrom = chrom = str2chrom(sx) ;

    if ((chrom>=MAXCH)  || (chrom <=0)) {
     if (nbad<10) printf("warning: bad chrom: %s", line) ;
     ++nbad ; 

     sdpt -> chrom = MIN(chrom, BADCHROM) ;  
     sdpt -> chrom = MAX(chrom, 0) ;  
     
     sdpt -> ignore = YES ;  
    }


    sdpt -> gpos = atof(spt[2]) ;
    if (sdpt->gpos > 100) {
      if (sdpt->gpos > 1.0e6) 
        fatalx("absurd genetic distance:\n%s\n", line) ;
     if (!usecm) {
      printf("*** warning.  genetic distances are in cM not Morgans\n") ;
      printf("%s\n",line) ;
     }
     usecm = YES ;
    }
    maxgpos[chrom] = MAX(maxgpos[chrom], sdpt -> gpos) ;
    maxg = MAX(maxg, maxgpos[chrom]) ;
    sdpt -> ppos = atof(spt[3]) ;
    if (nsplit<8) { 
     ivzero(sdpt->nn,4) ; 
    }
    else {
     for (k=0; k<4; k++) { 
      sdpt->nn[k] = atoi(spt[4+k]) ;
     }
    }
    sdpt -> inputrow = num ;
    ++num ;
   }
   freeup(spt, nsplit) ;
   continue ;
 }
 if (maxg<=0.00001)  {
  printf("genetic distance set from physical distance\n") ;
  usecm = NO ;
  for (k=0; k<num ; ++k) {
   snpraw[k] -> gpos = 1.0e-8 * snpraw[k] -> ppos ;
  }
 }
 if (usecm) {
  for (k=0; k<num ; ++k) {
   snpraw[k] -> gpos /= 100.0 ;
  }
 }
 if (pedcols == NULL) {
  ZALLOC(pedcols, num, int) ;  
  ivclear(pedcols, -1, num) ;
  numpedcols = num ;
 }
 
 fclose(fff) ;
 return num ;
}


int
numfakes(SNPDATA **snpraw, int *snpindx, int nreal, double spacing) 

// it seems better for this internal routine 
// to use the precomputed values
{
  int nignore, numsnps;
  int nfake = 0, i, k, indx ;
  int num=0; 
  SNP *cupt ;
  SNPDATA *sdpt ;
  char *sname ;
  int *sp ;
  int xc = 0, chrom ;
  double fakedis, realdis ;   // gpos for fake marker 
  double physpos ;

  if (spacing <= 0.0) fakedis = 1.0e20 ;

  for (k=0; k< nreal ; k++) {  

   indx = snpindx[k] ;
   sdpt = snpraw[indx] ;

   chrom = sdpt -> chrom ;
   realdis = sdpt -> gpos ;           

   if (chrom != xc)  {  
    fakedis = nextmesh(realdis, spacing) ;
    xc = chrom ;
   }
   while (fakedis<realdis) {  
     fakedis += spacing ;
     ++nfake ;
   }
  }
   return nfake ;
}
double nextmesh(double val, double spacing)  
{
 if (spacing==0.0) return 1.0e8 ;
 return ceil(val/spacing)*spacing ;
}


int
loadsnps(SNP **snpm, SNPDATA **snpraw, 
  int *snpindx, int nreal, double spacing, int *numignore) 

// snppos, snpindx could be recalculated but 
// it seems better for this internal routine 
// to use the precomputed values
// do NOT call externally
{
  int nignore, numsnps;
  int nfake = 0, i, k, indx ;
  int num=0; 
  SNP *cupt, *lastcupt ;
  SNPDATA *sdpt ;
  char *sname ;
  int *sp ;
  int xc = 0, chrom ;
  double fakedis, realdis, xrealdis ;   // gpos for fake marker 
  double physpos ;
  double xl, xr, xmid, al, ar, fraw ;
  int nn[2], n0, n1 ;
  int cnum ;
  int inputrow ;


  if (spacing <= 0.0) fakedis = 1.0e20 ;
  for (k=0; k< nreal ; k++) {  

   indx = snpindx[k] ;
   sdpt = snpraw[indx] ;

   chrom = sdpt -> chrom ;
   sname = sdpt -> ID ;
   realdis = sdpt -> gpos ;           
   physpos = sdpt -> ppos ;
   inputrow = sdpt -> inputrow ;

   if (chrom != xc)  {  
    fakedis = nextmesh(realdis, spacing) ;
    xc = chrom ;
    cnum = 0 ;
   }


   while (fakedis<realdis) {  

     if (cnum==0) break ;
     if (sdpt -> ignore) break ;

     if (nfake>=tempfake)  fatalx(" too many fake markers (bug) %d %d\n", num, nfake) ;
     if (num>=tempnum)  fatalx(" too many markers (bug) %d %d\n", num, nfake) ;
   
     cupt = snpm[num] ;
     if (cupt == NULL) fatalx("bad loadsnps\n") ;
     sprintf(cupt -> ID, "fake-%d:%d", xc, nfake) ;
     cupt -> estgenpos = cupt -> genpos = fakedis ; 
     xl = lastcupt -> genpos ;
     al =  lastcupt -> physpos ;
     xr =  realdis ;          ;
     ar =  physpos ;
     cupt -> physpos = interp(xl, xr, fakedis, al, ar) ;
     cupt -> markernum = num ;
     cupt -> isfake = YES ;
     cupt -> chrom  = xc ;
     fakedis += spacing ;
     ++num ;
     ++nfake ;
   }

   ++cnum ;
   cupt = snpm[num] ;
   if (cupt == NULL) fatalx("bad loadsnps\n") ;
   strcpy(cupt -> ID, sname) ;
   sdpt -> cuptnum = num ;
   cupt -> estgenpos = cupt -> genpos = realdis ;                          
   cupt -> physpos = physpos ;
   cupt -> markernum = num ;
   cupt -> isfake = NO ;
   cupt -> ignore  = sdpt -> ignore ;
   if ((cupt -> ignore == NO) && (cupt -> isfake == NO)) lastcupt = cupt ;
   cupt -> isrfake = sdpt -> isrfake ; 
   cupt -> chrom  = xc ;
   cupt -> tagnumber = inputrow ; // just used for pedfile 
   if (inputrow >=0) {
    if (inputrow >= numpedcols) fatalx("pedcols overflow\n") ;
    pedcols[inputrow] = num ;
   }

   n0 = sdpt->nn[0] ;
   n1 = sdpt->nn[1] ;
   fraw = mknn(nn, n0, n1) ;
   copyiarr(nn, cupt->af_nn, 2) ;
   cupt -> aftrue = cupt -> af_freq = fraw ;
   cupt -> aa_aftrue = cupt -> aa_af_freq = fraw ;

   if (sdpt -> alleles != NULL) { 
    cupt -> alleles[0] = sdpt -> alleles[0] ;
    cupt -> alleles[1] = sdpt -> alleles[1] ;
   }
   else {  
    cupt -> alleles[0] = '1' ; 
    cupt -> alleles[1] = '2' ; 
   }

   n0 = sdpt->nn[2] ;
   n1 = sdpt->nn[3] ;
   fraw = mknn(nn, n0, n1) ;
   copyiarr(nn, cupt->cauc_nn, 2) ;
   cupt -> cftrue = cupt -> cauc_freq = fraw ;
   cupt -> aa_cftrue = cupt -> aa_cauc_freq = fraw ;
   ++num ;
  }
// now make list of ignored snps used by loadgeno for check
  numsnps = num ;
  for (k=0; k< nreal ; k++) {  
   indx = snpindx[k] ;
   sdpt = snpraw[indx] ;
   if (sdpt->ignore == NO) continue ;
   inputrow = sdpt -> inputrow ;
   chrom = sdpt -> chrom ;
   sname = sdpt -> ID ;
   realdis = sdpt -> gpos ;           
   physpos = sdpt -> ppos ;
   cupt = snpm[sdpt -> cuptnum] ;
   cupt -> tagnumber = inputrow ; // just used for pedfile 
/**
   strncpy(cupt -> ID, sname, IDSIZE-1) ; 
   cupt -> genpos = realdis ;                          
   cupt -> physpos = physpos ;
   cupt -> markernum = num ;
   cupt -> isfake = NO ;
   cupt -> ignore = YES ;
   cupt -> chrom  = chrom ;
*/
   ++num ;
  }
  nignore = 0 ;
  for (k=0; k<numsnps; ++k) {
   cupt = snpm[k] ; 
   if (cupt -> ignore) ++nignore ;
  }
  *numignore = nignore ;
  return numsnps ;
}

double interp (double l, double r, double x, double al, double ar)
// linearly interp ;
{
  double y, y1, y2 ;
  y = (r-l) ;
  if (y==0.0) return 0.5*(al+ar) ;
  y1 = (r-x)/y ;
  y2 = (x-l)/y ;
  return y1*al +y2*ar ;
}

int getindivs(char *indivfname, Indiv ***indmarkpt) 
// returns number of indivs  
{

 static Indiv **indivmarkers ;
 int nindiv, i ;

 if (indivfname == NULL) fatalx("(getindivs) NULL indivfname\n") ;
 nindiv = getsizex(indivfname) ;
 if (nindiv <= 0) fatalx("no indivs found: indivname: %s\n", indivfname) ;
 ZALLOC(indivmarkers, nindiv, Indiv *) ;

 for (i=0; i<nindiv; i++) {  
  ZALLOC(indivmarkers[i], 1, Indiv) ;
 }
 clearind(indivmarkers, nindiv) ;
 *indmarkpt = indivmarkers ;
 readinddata(indivmarkers, indivfname) ;
 return nindiv ;
}

int readinddata(Indiv **indivmarkers, char *fname) 
{
 char line[MAXSTR] ;
 char *spt[MAXFF], *sx ;
 int nsplit, num=0, k ;
 int skipit ;
 Indiv *indx ;

 FILE *fff ;

 if (ispedfile(fname)) {  
  plinkinputmode = YES ;
  return  readindpeddata(indivmarkers, fname) ;
 }

 openit(fname, &fff, "r") ;
 while (fgets(line, MAXSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit == 0) continue ;
   sx = spt[0] ;
   skipit = setskipit(sx) ;
   if (skipit == NO) { 
    if (nsplit<3) fatalx("%s bad line: %s", fname, line) ;
    indx = indivmarkers[num] ;
    if (strlen(sx)>=IDSIZE) fatalx("ID too long: %s\n", sx) ;
    strcpy(indx->ID, sx) ;
    indx -> idnum = num ;
    sx = spt[1] ;
    indx -> gender = sx[0] ;
    indx -> affstatus = indx -> ignore = NO ;
    sx = spt[2] ;
    if (strcmp(sx, "Ignore") == 0) indx->ignore = YES ;
    if ((qtmode) && (!indx->ignore)) { 
     indx -> egroup = strdup("Case") ;
     indx -> qval = indx -> rawqval = atof(sx) ;
    }
    else {
     indx -> egroup = strdup(sx) ;  
    }
// affstatus set by setstatus  
    ++num ;
   }
   freeup(spt, nsplit) ;
   continue ;
 }

 fclose(fff) ;
 return num ;
}

int readindpeddata(Indiv **indivmarkers, char *fname) 
///
{
 char *line ;
 char *spt[MAXFF], *sx, *sx0, gender ;
 int nsplit, num=0, k, i ;
 int skipit ;
 Indiv *indx ;
 int nindiv ; 
 int maxnsplit = 0 ;
 char nnbuff[IDSIZE] ;
 int nok = 0 ;

 FILE *fff ;

 maxgenolinelength = maxlinelength(fname) ;
 ZALLOC(line, maxgenolinelength+1, char) ;
 openit(fname, &fff, "r") ;

 while (fgets(line, maxgenolinelength, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit == 0) continue ;
   sx0 = sx = spt[0] ;
   skipit = NO ;
   if (sx[0] == '#') skipit = YES ;
   if (skipit == NO) { 
    if (nsplit<6) fatalx("%s bad line: %s", fname, line) ;
    indx = indivmarkers[num] ;
    if (strlen(sx)>=IDSIZE) fatalx("ID too long: %s\n", sx) ;
    maxnsplit = MAX(maxnsplit, nsplit) ;

    sx = spt[1] ;
    pedname(nnbuff, sx0, sx) ;
    strcpy(indx->ID, nnbuff) ;
    indx -> idnum = num ;
    sx = spt[4] ;
    k = atoi(sx) ;
    gender = 'U' ;
    if (k==1) gender = 'M' ;
    if (k==2) gender = 'F' ;
    indx -> gender = gender ;
    indx -> affstatus = indx -> ignore = NO ;

    sx = spt[5] ;
    if (qtmode) {
     indx -> egroup = strdup("Case") ;
     indx -> qval = indx -> rawqval = atof(sx) ;
    }

    else {
     k = 99 ;
     if (strcmp(sx, "-9") == 0) k = -9 ;
     if (strcmp(sx, "9") == 0) k = 9 ;
     if (strcmp(sx, "0") == 0) k = 0 ;
     if ((pedignore == NO) && (k==0)) k = 3 ;
     if (strcmp(sx, "1") == 0) k = 1 ;
     if (strcmp(sx, "2") == 0) k = 2 ;
     switch (k) { 
      case 9:  
        indx -> ignore = YES ;
        printf("%s ignored\n", indx -> ID) ;
        break ;
      case -9:  
        indx -> ignore = YES ;
        printf("%s ignored\n", indx -> ID) ;
        break ;
      case 0:  
        indx -> ignore = YES ;
        printf("%s ignored\n", indx -> ID) ;
        break ;
      case 1: 
       indx -> egroup = strdup("Control") ;
       break ; 
      case 2: 
       indx -> egroup = strdup("Case") ;
       break ;
      case 3: 
       indx -> egroup = strdup("???") ;
       break ;
      default:  
       indx -> egroup = strdup(sx) ;
     }

    }

// affstatus set by setstatus  
    if (indx -> ignore == NO) ++nok ;
    ++num ;
   }
   freeup(spt, nsplit) ;
   continue ;
 }

 if (nok == 0) {
   printf("all individuals set ignore.  Likely input problem (col 6)\n") ;
   printf("resetting all individual...\n") ;
   for (i=0; i<num; i++) {  
    indx = indivmarkers[i] ;
    indx -> ignore = NO ;
    indx -> egroup = strdup("???") ;
   }
  }
 
   

 if (maxnsplit<8) maxgenolinelength = -1 ;
 free(line) ;

 fclose(fff) ;
 return num ;
}
void pedname(char *cbuff, char *sx0, char *sx1) 
{
  int l0 , l1, ll ;

  l0 = strlen(sx0) ;
  l1 = strlen(sx1) ;
  ll = l0 +l1 + 1 ; 
  if (familynames == NO) ll = l1 ;
  if (ll>=IDSIZE) {
   fatalx("idnames too long %s %s ll: %d limit: %d\n", sx0, sx1, ll, IDSIZE-1) ;
  }
  if (familynames == YES) {
   strcpy(cbuff, sx0) ; 
   cbuff[l0] = ':' ; 
   strcpy(cbuff+l0+1, sx1) ;
   return ;
  }
  strcpy(cbuff, sx1) ;

}
  
int readtldata(Indiv **indivmarkers, int numindivs, char *inddataname)  
{ 
 char line[MAXSTR] ; 
 char *spt[MAXFF], *sx ; 
 int nsplit, num=0, k, ind ; 
 int skipit ; 
 Indiv *indx ; 
 double y ; 
 double gg[3] ;  
 
 FILE *fff ; 
 openit(inddataname, &fff, "r") ; 
 while (fgets(line, MAXSTR, fff) != NULL) { 
   nsplit = splitup(line, spt, MAXFF) ;  
   if (nsplit == 0) continue ; 
   sx = spt[0] ; 
   skipit = NO ; 
   if (strcmp(sx, "Indiv_Index") == 0) {
// hack.  thetafile should be output with leading ## 
    freeup(spt, nsplit) ; 
    continue ; 
   }
   if (sx[0] == '#') {                
    freeup(spt, nsplit) ; 
    continue ; 
   }
    if (nsplit<8) fatalx("%s bad line: %s", inddataname, line) ; 
    sx = spt[1] ;
    ind = indindex(indivmarkers, numindivs, sx)  ; 
    if (ind<0) fatalx("(readtldata) indiv: %s not found \n", sx) ; 
    indx = indivmarkers[ind] ; 

    indx -> theta_mode = atof(spt[3]) ;
    indx -> lambda_mode  = atof(spt[7]) ;
    indx -> Xtheta_mode = atof(spt[5]) ;
    indx -> Xlambda_mode  = atof(spt[9]) ;

    freeup(spt, nsplit) ; 
    continue ; 
 } 
  
 fclose(fff) ; 
 return num ; 
} 

int readfreqdata(SNP **snpm, int numsnps, char *inddataname)  
{ 
 char line[MAXSTR] ; 
 char *spt[MAXFF], *sx ; 
 int nsplit, num=0, k, ind ; 
 int skipit ; 
 SNP *cupt ; 
 
 FILE *fff ; 
 openit(inddataname, &fff, "r") ; 
 while (fgets(line, MAXSTR, fff) != NULL) { 
   nsplit = splitup(line, spt, MAXFF) ;  
   if (nsplit == 0) continue ; 
   sx = spt[0] ; 
   skipit = NO ; 
   if (sx[0] == '#') {                
    freeup(spt, nsplit) ; 
    continue ; 
   }
    if (nsplit<6) fatalx("%s bad line: %s", inddataname, line) ; 
    sx = spt[2] ;
    ind = snpindex(snpm, numsnps, sx)  ; 
    if (ind<0) fatalx("(readfreqdata) snp %s not found \n", sx) ; 
    cupt = snpm[ind] ; 
    cupt -> aa_af_freq = cupt -> af_freq  = atof(spt[3]) ;
    cupt -> aa_cauc_freq = cupt -> cauc_freq  = atof(spt[5]) ;

    freeup(spt, nsplit) ; 
    ++num ;
    continue ; 
 } 
  
 fclose(fff) ; 
 return num ; 
} 



int setstatus(Indiv **indm, int numindivs, char *smatch)   
// return number set 
// smatch = NULL => set everything 
{ 
  return setstatusv(indm, numindivs, smatch, YES) ;
}

int setstatusv(Indiv **indm, int numindivs, char *smatch, int val)   
// return number set 
// smatch = NULL => set everything 
{ 
  int i, n=0 ; 
  Indiv *indx ; 
  char *sx ; 
  for (i=0; i<numindivs; i++) { 
   indx = indm[i] ; 
   if (indx -> ignore) continue ;
   sx = indx -> egroup ; 
   if (smatch == NULL) {
    ++n ; 
    indx->affstatus = val ; 
    continue ;
   } 
   if (strcmp(sx, smatch) == 0) {  
    ++n ; 
    indx->affstatus += val ; 
   } 
  } 
  return n ; 
} 


long getgenos(char *genoname, SNP **snpmarkers, Indiv **indivmarkers, 
 int numsnps, int numindivs, int nignore) 
{
// read genofile.  Use hashtable to improve search 
// if genofile is gzipped decompress to trashdir
 char *gname, *genotmp = NULL ;
 ENTRY *hashlist, *iteml  ;
 ENTRY item1 ;
 int k, num, indiv, lgt ;
 int val ;
 void *basept = 0 ;
 int bigoff ;
 int tcheck ;

// we use a trick: want to store k 
// store basept + k instead 
// basept + k + bigoff for individual ID 

 SNP *cupt ;
 Indiv *indx ;

 char line[MAXSTR], cmd[MAXSTR] ;
 char *spt[MAXFF], *sx ;
 int nsplit, nsnp ;
 int skipit, kret, tpackmode, teigenstratmode ;

 FILE *fff ;
 int gnlen, ngenos=0 ;

 double y ;
 char *pbuff ;

 item1.key  = NULL ;           
 item1.data = NULL ;           

 if (genoname == NULL) fatalx("(getgenos) NULL genoname\n") ; 
 gname = genoname ;
 gnlen = strlen(genoname) ;
 if (strcmp(genoname+gnlen-3, ".gz") == 0)  {  
  makedir(trashdir) ; 
  sprintf(line, "%s/genotmp:%d", trashdir, getpid()) ;
  genotmp = strdup(line) ;
  sprintf(cmd, "gunzip -c %s > %s", genoname, genotmp) ;
  printf("unzip cmd: %s\n", cmd) ;
  system (cmd) ;
  kret = system (cmd) ;
  if (kret<0) {  
   perror("gunzip failed\n") ;
   fatalx("gunzip failed... probably out of disk space\n") ;
  }
  printf("geno file unzipped\n") ;
  gname = genotmp ;
 }

 tcheck = checksize(numsnps, numindivs, outputmode)  ;     
 if (tcheck == -2) fatalx("Data sets with more than 8 billion genotypes are not permitted\n") ;
 if (tcheck == -1) fatalx("Output files of size >2GB are not permitted: use a more compact output data format. Also see documentation of chrom, badsnpname and checksizemode parameters.\n") ;

 if (ispedfile(gname))   {

   if (pedcols == NULL) fatalx("pedcols not allocated (no map file ?)") ;
   getpedgenos(genoname, snpmarkers, indivmarkers, 
     numsnps,  numindivs,  nignore)  ;
   freeped()  ;
   return ; 
 }  

 if (isbedfile(gname))   {
   return  getbedgenos(genoname, snpmarkers, indivmarkers, 
     numsnps,  numindivs,  nignore)  ;
 }  

 tpackmode = ispack(gname) ;
 

 nsnp = numsnps  ;
 if (tpackmode) { 
   inpack(gname, snpmarkers, indivmarkers,  nsnp, numindivs)  ; 
   for (k=0; k<nsnp; k++)  { 
    cupt = snpmarkers[k] ;
    if (cupt -> ignore) continue ;
    if ((cupt -> isfake) && (!(cupt -> isrfake))) continue ;
    cupt -> ngtypes = numindivs ; 
    if (cupt -> gtypes == NULL) ZALLOC(cupt -> gtypes, 1, int) ; 
   }
   packmode = YES ;
   return nsnp*numindivs ;
 }
 teigenstratmode = iseigenstrat(gname) ;
 if (teigenstratmode) { 
  packmode = YES ;
  ineigenstrat(gname, snpmarkers, indivmarkers,  nsnp, numindivs)  ; 
  freeped() ;
  return nsnp*numindivs ;
 }

 

 y = (double) (numindivs * 2) / (8 * (double) sizeof (char)) ;  
 rlen = nnint(ceil(y)) ;
 packlen = rlen*numsnps ;
 if (packmode) {
   ZALLOC(packgenos, packlen, char) ;
   pbuff = packgenos ;
   clearepath(packgenos) ;
 }

 num = nsnp + numindivs ;
 xhcreate(5*num) ;
 ZALLOC(hashlist, num, ENTRY) ;
 bigoff = nsnp + 100 ;

 for (k=0; k<nsnp; k++)  { 
  cupt = snpmarkers[k] ;
  if ((cupt -> isfake) && (!(cupt -> isrfake))) continue ;
  iteml = hashlist+k ;
  iteml -> key = cupt->ID  ;
  iteml -> data = basept +k ;  
  if (xhsearch(*iteml, FIND)  != NULL) 
   fatalx("duplicate ID: %s\n", iteml -> key) ;
  (void) xhsearch(*iteml, ENTER) ;
 }

 for (k=0; k<numindivs; k++)  { 
  indx = indivmarkers[k] ;
  iteml = hashlist+numsnps+k ;
  iteml -> key = indx->ID  ;
  iteml -> data = basept+k+bigoff ;  
  if (xhsearch(*iteml, FIND)  != NULL) 
   fatalx("duplicate ID: %s\n", iteml -> key) ;
  (void) xhsearch(*iteml, ENTER) ;
 }
 openit(gname, &fff, "r") ;
 while (fgets(line, MAXSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit == 0) continue ;
   sx = spt[0] ;
   skipit = NO ;
   if (sx[0] == '#') skipit = YES ;
   skipit = setskipit(sx) ;
   if (skipit == NO) { 
    if (nsplit<3)  fatalx("bad geno line.  missing field?\n", line) ;
    item1.key = spt[0] ;
    iteml = xhsearch(item1, FIND) ; 
    if (iteml == NULL) {
     fatalx("(genotypes) bad ID (SNP): %s\n", line) ;
    }
    k = (int) (iteml->data -  basept) ;

    if (k>=numsnps) {
      fatalx("bad genotype line: `snp' may be Indiv Id\n%s\n", line) ;
    }

    cupt = snpmarkers[k] ;
    if (cupt -> ignore) {
      freeup(spt, nsplit) ;
      continue ;
    }
    item1.key = spt[1] ;
    iteml = xhsearch(item1, FIND) ; 
    if (iteml == NULL) {
     fatalx("(genotypes) bad ID: (Indiv) %s\n", line) ;
    }
    indiv = (int) (iteml->data -  basept) ;
    indiv -= bigoff ;
    val = atoi(spt[2]) ;

    indx = indivmarkers[indiv] ;
    if (indx->ignore)  val = -1 ;
    if (checkxval(cupt, indx, val) == NO) val = -1 ;
    if (val>2) {  
     printf("*** warning invalid genotype: %s %s %d\n", 
       cupt -> ID, indx -> ID, val) ;  
       val = -1 ;
    }

    if (cupt -> ngtypes == 0) { 
     if (packmode == NO) {
      ZALLOC(cupt -> gtypes, numindivs, int) ;
     }
     else  {
      ZALLOC(cupt -> gtypes, 1, int) ;
      cupt -> pbuff = pbuff ;  
      pbuff += rlen ;
     }
     cupt -> ngtypes = numindivs ;
     for (k=0; k<numindivs; ++k) {
      putgtypes(cupt, k, -1) ;
     }
    }
    putgtypes(cupt, indiv, val) ;
    ++ngenos ;
   }
   freeup(spt, nsplit) ;
 }
 fclose(fff) ;
 free(hashlist) ;
 xhdestroy() ;
 if (genotmp != NULL) { 
   unlink(gname) ;
 }
 /* printf("genotype file processed\n") ; */
 freeped() ; 
 return ngenos ;
}

void freeped() 
// destructor for pedcols 
{
 if (pedcols == NULL) return ; 
 if (dofreeped == NO) return ;
 free(pedcols) ;
 pedcols = NULL ;
 numpedcols = 0 ;
 maxgenolinelength = -1 ;
}

int checkxval(SNP *cupt, Indiv *indx, int val) 
// check Male X marker not het
{
   if (cupt -> chrom != 23) return YES ;
   if (indx -> gender != 'M') return YES ;
   if (val != 1) return YES ;  
   if (malexhet)  return YES ;
   return NO; 
}

void clearsnp(SNP *cupt) 
{

   cupt -> af_freq =
   cupt -> cauc_freq = -1 ;
   cupt -> aa_af_freq = 
   cupt -> aa_cauc_freq = -1 ; 
   cupt -> estgenpos = 0 ;
   cupt -> genpos = 0 ;
   cupt -> physpos = 0 ;
   cupt -> ngtypes = 0 ;
   cupt -> pbuff = NULL ;
   cupt -> ebuff = NULL ;
   cupt -> gtypes = NULL ;
   cupt -> modelscores = NULL ;
   cupt -> totmodelscores = NULL ;
   cupt -> score = cupt -> weight = 0.0 ;
   cupt -> isfake = NO ; 
   cupt -> ignore = NO ; 
   cupt -> isrfake = NO ; 
   cupt -> estdis = 0 ; 
   cupt -> dis = 0 ; 
   cupt -> esum = 0 ; 
   cupt -> lsum = 0 ;
   cupt -> gpsum = 0 ; 
   cupt -> gpnum = 0 ;
   cupt -> pcupt = NULL ;
   cupt -> tagnumber = -1 ;
   cclear((unsigned char *) cupt -> alleles, CNULL, 2)  ;
}

int rmindivs(SNP **snpm, int numsnps, Indiv **indivmarkers, int numindivs) 
// squeeze out ignore
// dangerous bend.  Of course indivmarkers indexing will change
{
   int n = 0, g, i, k ;
   int x ; 
   Indiv *indx ;
   SNP *cupt ;

   for (k=0; k<numindivs; ++k) {  
    if (indivmarkers[k] -> ignore == YES)  continue ;
    if (n==k) { 
     ++n ;
     continue ; 
    }
// copy k -> n 
    indx = indivmarkers[n] ;
    indivmarkers[n] = indivmarkers[k] ;
    indx -> idnum = n ;
    for (i=0; i<numsnps; i++) {  
     cupt = snpm[i] ;
     if (cupt -> ignore) continue ;
     g = getgtypes(cupt, k) ;  
     putgtypes(cupt, n, g) ;
    }
    ++n ;
   }

   for (i=0; i<numsnps; i++) {  
    cupt = snpm[i] ;
    cupt -> ngtypes = n ; 
   }
   return n ;

}

int rmsnps(SNP **snpm, int numsnps) 
{
   int i,x ;
   SNP *cupt ;
   int lastc, chrom ;

// wipe out fakes not between real markers
   freesnpindex() ;
   lastc = -1 ;
   for (i=0; i<numsnps; i++)  {  
     cupt = snpm[i] ;
     if (cupt -> ignore) continue ;  
     chrom = cupt -> chrom ;  
     if ( (cupt -> isfake) && (chrom != lastc)) cupt -> ignore = YES ;
     if (!cupt -> isfake) lastc = chrom ;
   }

   lastc = -1 ;
   for (i=numsnps-1; i>=0; i--)  {  
     cupt = snpm[i] ;
     if (cupt -> ignore) continue ;  
     chrom = cupt -> chrom ;  
     if ( (cupt -> isfake) && (chrom != lastc)) cupt -> ignore = YES ;
     if (!cupt -> isfake) lastc = chrom ;
   }



   x = 0 ;
   for (i=0; i<numsnps; i++)  {  
    cupt = snpm[i] ;
    if (cupt -> ignore) {
     freecupt(&cupt) ;
     continue ; 
    }
    snpm[x] = snpm[i] ;
    ++x ;
   }

   for (i=0; i<x; i++)  {  
    cupt = snpm[i] ;
    cupt -> markernum = i ;
   }

   return x ;
}
void freecupt(SNP **cuppt) 
{
   SNP *cupt ; 
   cupt = *cuppt ;
    if (cupt -> modelscores != NULL) { 
     free(cupt->modelscores) ;
    }
    if (cupt -> totmodelscores != NULL) { 
     free(cupt->totmodelscores) ;
    }
    free(cupt) ;
    cupt = NULL ;
}

void
clearind(Indiv **indm, int numind)   
{
   Indiv *indx ;  
   double theta ; 
   int i ;

   for (i=0; i<numind; i++)  {
    indx = indm[i] ;
    indx -> egroup = NULL ;
    indx -> affstatus = indx -> ignore = NO ;
    indx -> gender = 'U' ;
    indx = indm[i] ;
    indx -> Xtheta_mode = indx->theta_mode = a1/(a1+b1) ;
    indx -> Xlambda_mode = indx -> lambda_mode = lp1/lp2 ;
    indx -> thetatrue = -1.0 ;  // silly value
    indx -> qval = indx -> rawqval = 0.0 ;
   }
   cleartg(indm, numind) ;
}

void cleartg(Indiv **indm, int nind) 
{
    int i ;
    Indiv *indx ;

      for (i=0; i< nind; i++) {  
       indx = indm[i] ;
       vzero(indx -> totgamms, 3) ;
       indx -> totscore = 0.0 ; 
      }

}


double mknn(int *nn, int  n0, int n1) 
{

  double x ;
  int t ;

  nn[0] = n0 + 1 ;       
  nn[1] = n1 + 1 ;       

// no clipping.  (Old code clipped here)
  t = intsum(nn,2) ;
  x = ((double) nn[0]) / (double) t ;

  return x ;
}

void
setug(Indiv **indm, int numind, char gender)   
{

   Indiv *indx ;  
   double theta ; 
   int i ;

   for (i=0; i<numind; i++)  {
    indx = indm[i] ;
    if (indx -> gender == 'U') indx -> gender = gender ;
   }

}


void dobadsnps(SNPDATA **snpraw, int nreal, char *badsnpname) 
{

  FILE *fff ;
  char line[MAXSTR] ;
  char *spt[MAXFF] ;
  char *ss ;
  int indx, nsplit, n ;

  if (badsnpname == NULL) return ;
  openit (badsnpname, &fff, "r") ;

  while (fgets(line, MAXSTR, fff) != NULL)  {
    nsplit = splitup(line, spt, MAXFF) ; 
    if (nsplit==0) continue ;
    if (spt[0][0] == '#') { 
     freeup(spt, nsplit) ;
     continue ;
    }
    indx = snprawindex(snpraw, nreal, spt[0]) ;
    if (indx >=0) {
     snpraw[indx] -> ignore = YES ;
     if ((nsplit >=2) && (checkfake(spt[1]))) {  
       snpraw[indx] -> ignore = NO ;
       snpraw[indx] -> isrfake = YES ;
     }
    }
    freeup(spt, nsplit) ;
  }
  fclose (fff) ;
}
int checkfake(char *ss) 
// yes if string ss is "Fake" 
// ss is overwritten
{
  ss[0] = tolower(ss[0]) ;
  if (strcmp(ss,"fake") == 0) return YES ;
  return NO ;

}

void printsnps(char *snpoutfilename, SNP **snpm, int num, 
  Indiv **indm, int printfake, int printvalids) 
{ 
  int i, ppos ;
  SNP *cupt  ;
  char ss[MAXSTR] ;
  FILE *xfile ;
  int numvcase, numvcontrol ;
  char c ;

  if ((snpoutfilename != NULL) && (strcmp(snpoutfilename, "NULL") == 0))  return ;
  if (snpoutfilename != NULL)  {
   openit(snpoutfilename, &xfile, "w") ;
  }
  else xfile = stdout ;

  if( tersemode == NO) {
    fprintf(xfile,"\n");
    fprintf(xfile,"###DETAILS ABOUT THE MARKERS\n");
    fprintf(xfile,"##Gen_Pos: genetic position, Phys_pos: Physical position\n");
    fprintf(xfile,"##Afr_vart: Parental African variant allele count, Afr_ref: Parental African reference allele count\n");
    fprintf(xfile,"##Eur_vart: Parental European variant allele count, Eur:ref:Parental European reference allele count\n");   

    fprintf(xfile, "\n") ;
    fprintf(xfile,"%20s %5s %10s %18s",  "#SNP_Id","Chr_Num","Gen_Pos","Phys_Pos") ;
    fprintf(xfile,  " %9s %9s %9s %9s" ,"Afr_vart","Afr_ref","Eur_vart","Eur_ref");
    fprintf(xfile, "\n") ;
  }
  for (i = 0; i <num; ++i) {  
     cupt = snpm[i] ;
     if (outputall==NO)  {
      if (!printfake && (ignoresnp(cupt))) continue ;
      if (!printfake && (cupt -> isrfake)) continue ;
     }
     fprintf(xfile, "%20s %3d ", cupt->ID, cupt -> chrom) ;
     fprintf(xfile, "%15.6f %15.0f",  cupt -> genpos, cupt -> physpos) ;

     if (tersemode) { 
      printalleles(cupt, xfile) ;
      fprintf(xfile, "\n") ;
      continue ;
     }

     fprintf(xfile, " %8d ", cupt -> af_nn[0]) ;
     fprintf(xfile, "%8d ", cupt -> af_nn[1]) ;
     fprintf(xfile, "%8d ", cupt -> cauc_nn[0]) ;
     fprintf(xfile, "%8d", cupt -> cauc_nn[1]) ;
     if (!printvalids) { 
      printalleles(cupt, xfile) ;
      fprintf(xfile, "\n") ;
      continue ;
     }
     numvcase = numvalidgtx(indm, cupt, 1) ;
     numvcontrol = numvalidgtx(indm, cupt, 0)  ;
     fprintf(xfile, "   %6d %6d",numvcase, numvcontrol) ;
     fprintf(xfile, "    %d %d %d", cupt -> ignore, cupt -> isfake, cupt -> isrfake) ;
     printalleles(cupt, xfile) ;
     fprintf(xfile, "\n") ;
  }
  if (snpoutfilename != NULL) 
    fclose(xfile) ;
}
void printalleles(SNP *cupt, FILE *fff)  
{

  char c ;
  if ((c = cupt -> alleles[0]) != CNULL) fprintf(fff, " %c", c) ;
  if ((c = cupt -> alleles[1]) != CNULL) fprintf(fff, " %c", c) ;

}

void printdata(char *genooutfilename, char *indoutfilename,
  SNP **snpm, Indiv **indiv, int numsnps, int numind, int packem)
{

  FILE *gfile, *ifile;
  int i,j, t;
  SNP *cupt;
  Indiv *indx;
  char ss[MAXSTR];
  char *gfilename ;
  int dogenos = YES ;

  if (packem) 
   printf("packedancestrymap output\n") ;
  else 
   printf("ancestrymap output\n") ;

  if ((genooutfilename != NULL) && (strcmp(genooutfilename, "NULL") == 0))  dogenos = NO ;
  if (genooutfilename == NULL) dogenos = NO ;

  if (dogenos) {
   gfilename = genooutfilename ;
   if (packem) { 
    outpack(genooutfilename, snpm, indiv, numsnps, numind) ;
     gfilename = NULL ;
    }

   if (gfilename != NULL)  {
    openit(gfilename, &gfile, "w") ;
    if(tersemode == NO) fprintf(gfile,"#SNP_ID,INDIV_ID,VART_ALLELE_CNT\n");
   }

   for (i = 0; i< numsnps; i++)  {
      if (gfilename == NULL) break ;
      cupt= snpm[i];

      if (outputall==NO)  {
       if (ignoresnp(cupt)) continue ;
       if (cupt -> isrfake) continue ;
      }

      for(j=0; j < cupt->ngtypes; j++)
	{
	  indx = indiv[j];
          if (indx -> ignore) continue  ;  
	  fprintf(gfile,"%20s %20s %3d\n",cupt->ID,indx->ID, getgtypes(cupt,j)) ;
	}
     }
   

   if(gfilename != NULL)
    fclose(gfile);

  }
  
  if(indoutfilename == NULL)
    return;
  if ((indoutfilename != NULL) && (strcmp(indoutfilename, "NULL") == 0))  return ;
  if (indoutfilename != NULL) 
    openit(indoutfilename, &ifile, "w") ;

  /* fprintf(ifile,"#INDIV,GENDER,POPULATION\n"); */
  for(i = 0; i< numind; i++)
    {
      indx = indiv[i];
      if (indx->ignore) continue ;            
      strcpy(ss, indx -> egroup) ;
      if ((qtmode) && (!indx->ignore)) {  
        sprintf(ss, "%9.3f", indx -> rawqval) ;
      }
      if (tersemode) { 
       fprintf(ifile,"%20s %c %10s",indx->ID, indx->gender,ss);
       fprintf(ifile,"\n") ;
       continue ;
      }
      t = numvalids(indx, snpm, 0, numsnps-1) ;
      fprintf(ifile,"%20s %c %10s %5d\n",indx->ID, indx->gender,ss, t);
    }

  if(indoutfilename != NULL)
    fclose(ifile);

}

int readindval(Indiv **indivmarkers, int numindivs, char *inddataname)  
{ 
 char line[MAXSTR] ; 
 char *spt[MAXFF], *sx ; 
 int nsplit, num=0, k, ind ; 
 int skipit ; 
 Indiv *indx ; 
 double y ; 
 double gg[3] ;  
 
 FILE *fff ; 
 openit(inddataname, &fff, "r") ; 
 for (k=0; k <numindivs; ++k) { 
    indx = indivmarkers[k] ; 
    indx -> affstatus = NO ;
    indx -> qval = -999.0 ;
 }
 while (fgets(line, MAXSTR, fff) != NULL) { 
   nsplit = splitup(line, spt, MAXFF) ;  
   if (nsplit < 2) {
    freeup(spt, nsplit) ; 
    continue ; 
   }
   sx = spt[0] ; 
   if (strcmp(sx, "Indiv_Index") == 0) {
    freeup(spt, nsplit) ; 
    continue ; 
   }
   if (sx[0] == '#') {                
    freeup(spt, nsplit) ; 
    continue ; 
   }
    ind = indindex(indivmarkers, numindivs, sx)  ; 
    if (ind<0) fatalx("(readindval) indiv: %s not found \n", sx) ; 
    indx = indivmarkers[ind] ; 
    indx -> qval = atof(spt[1]) ;
    indx -> affstatus = YES ;
    freeup(spt, nsplit) ; 
    continue ; 
 } 
  
 fclose(fff) ; 
 return num ; 
} 


int readgdata(Indiv **indivmarkers, int numindivs, char *gname)  
     // only needed for logreg 
     // not correct for X chromosome 
     // Needs correction  for males 
{ 
 char line[MAXSTR] ; 
 char *spt[MAXFF], *sx ; 
 int nsplit, num=0, k, ind ; 
 int skipit ; 
 Indiv *indx ; 
 double y ; 
 double gg[3] ;  
 
 FILE *fff ; 

 cleartg(indivmarkers, numindivs) ;
 openit(gname, &fff, "r") ; 
 while (fgets(line, MAXSTR, fff) != NULL) { 
   nsplit = splitup(line, spt, MAXFF) ;  
   if (nsplit == 0) continue ; 
   sx = spt[0] ; 
   skipit = NO ; 
   skipit = setskipit(sx) ;
   if (skipit == NO) {  
    if (nsplit<4) fatalx("%s bad line: %s", gname, line) ; 
    ind = indindex(indivmarkers, numindivs, sx)  ; 
    if (ind<0) fatalx("(readgdata) indiv: %s not found \n", sx) ; 
    indx = indivmarkers[ind] ; 
    for (k=0; k<3; k++) {  
     gg[k] = atof(spt[k+1]) ; 
    } 
    y = asum(gg, 3) ;  
    vst(gg, gg, 1.0/y, 3) ; 
    y = 0.5*(gg[1]+2.0*gg[2]) ;  /* est caucasian ancestry */ 
    indx -> thetatrue = y ;   
    copyarr(gg, indx -> totgamms, 3) ; 
   } 
   freeup(spt, nsplit) ; 
   continue ; 
 } 
  
 fclose(fff) ; 
 return num ; 
} 
int getweights(char *fname, SNP **snpm, int numsnps)  
// number of real lines 
{
 char line[MAXSTR] ;
 char *spt[MAXFF], *sx ;
 int nsplit, num=0 ;
 int skipit, k ;
 double weight ;

 FILE *fff ;
 for (k=0; k<numsnps; ++k) { 
  snpm[k] -> weight = 1.0 ;
 }
 openit(fname, &fff, "r") ;
 while (fgets(line, MAXSTR, fff) != NULL)  {
   nsplit = splitup(line, spt, MAXFF) ; 
   if (nsplit == 0) { 
     continue ;
   }
   sx = spt[0] ;
   skipit = NO ;
   skipit = setskipit(sx) ;
   k = snpindex(snpm, numsnps, sx) ;
   if (k<0) skipit = YES ;
   if (skipit == NO) {
    if (nsplit >1) {  
     sx = spt[0] ;
     weight =  atof(sx) ;
     snpm[k] -> weight = weight ;
     printf("weight set: %20s %9.3f\n", snpm[k] -> ID, weight) ;
     ++num ;
    }
   }
   freeup(spt, nsplit) ;
   continue ;
 }
 fclose(fff) ;
 fflush(stdout) ;
 return num ;
}
void
outpack(char *genooutfilename, SNP **snpm, Indiv **indiv, int numsnps, int numind) 
{
   char **arrx ;  
   int n, num, ihash, shash, i, g, j, k ;
   int nind , nsnp, irec ;
   Indiv *indx ; 
   SNP *cupt ;
   double y ;
   unsigned char *buff  ;
   int fdes, ret ;
   char *packit ;

   n = numind ;
   ZALLOC(arrx, n, char *) ;

   num = 0 ;
   for (i=0; i<n ; i++)  {  
    indx = indiv[i] ;  
    if ((outputall == NO ) && indx -> ignore) continue ;
    arrx[num] = strdup(indx -> ID) ;
    ++num ;
   }
   ihash = hasharr(arrx, num) ;
   nind= num ;


   freeup(arrx, num) ;  
   free(arrx) ;
   
   n = numsnps ;
   ZALLOC(arrx, n, char *) ;
   num = 0 ;
   for (i=0; i<n ; i++)  {  
    cupt = snpm[i] ;  
    if (outputall==NO)  {
     if (ignoresnp(cupt)) continue ;
     if (cupt -> isrfake) continue ;
    }
    arrx[num] = strdup(cupt -> ID) ;
    ++num ;
   }
   shash = hasharr(arrx, num) ;
   nsnp = num ;
   freeup(arrx, num) ;  
   free(arrx) ;

   // printf("ihash:  %x   shash: %x\n", ihash, shash) ; 
   y = (double) (nind * 2) / (8 * (double) sizeof (char)) ;  
   rlen = nnint(ceil(y)) ;
   rlen = MAX(rlen, 48)  ;  
   // printf("nind: %d  rlen:  %d\n", nind, rlen) ;
   ZALLOC(buff, rlen, unsigned char) ;
   sprintf((char *) buff,"GENO %7d %7d %x %x", nind, nsnp, ihash, shash) ;

   ridfile(genooutfilename) ;
   fdes = open(genooutfilename, O_CREAT | O_TRUNC | O_RDWR, 0666);
 
   if (fdes<0) {
    perror("bad genoout") ;
    fatalx("open failed for %s\n", genooutfilename) ;
   }
   if (verbose) 
    printf("file %s opened\n", genooutfilename) ;

    ret = write(fdes, buff, rlen) ;
    if (ret<0) {  
     perror("write failure") ;
     fatalx("(outpack) bad write") ;
    } 

    irec = 1; 
    for (i=0; i<numsnps ; i++)  {  
     cupt = snpm[i] ;  
      if (outputall==NO)  {
       if (ignoresnp(cupt)) continue ;
       if (cupt -> isrfake) continue ;
      }
     cclear((unsigned char *) buff, 0X00, rlen) ;  
     num = 0 ;
     for (j=0; j< numind; j++) { 
      indx = indiv[j] ;  
      if (indx -> ignore) continue ;
      g = getgtypes(cupt, j) ;
      if (g<0) g=3 ;
      wbuff( buff, num, g) ;
      ++num ;
     }
     ret = write(fdes, buff, rlen) ;
     if (ret<0) {  
      perror("write failure") ;
      fatalx("(outpack) bad write") ;
     } 
     if (verbose) {
      printf("record: %4d ", irec) ; 
      for (k=0; k<rlen; ++k) {  
       printf(" %02x", (unsigned char) buff[k]) ;
      }
      printf("\n") ;
     }
     ++irec ;
    }
    close(fdes) ;
    free(buff) ;
    // printf("check: %s %d\n", genooutfilename, ispack(genooutfilename)) ;
}
int ispack(char *gname) 
// checks if file is packed gfile 
{
  int fdes, t, ret ;  
  char buff[8] ;

  fdes = open(gname, O_RDONLY) ;
  if (fdes<0) {  
   perror("open failure") ;
   fatalx("(ispack) bad open %s\n", gname) ;
  } 
  t = read(fdes, buff, 8 ) ;
  if (t<0) {  
   perror("read failure") ;
   fatalx("(ispack) bad read") ;
  } 
  close(fdes) ;
  buff[4] = '\0' ;  
  ret = strcmp(buff, "GENO") ;
  if (ret == 0) return YES ;
  return NO ;
  
}
int iseigenstrat(char *gname)  

{

   FILE *fff ;
   char line[MAXSTR] ;
   char *spt[MAXFF], *sx ;
   int nsplit, num=0, k ;


   openit(gname, &fff, "r") ;

   while (fgets(line, MAXSTR, fff) != NULL)  {
    nsplit = splitup(line, spt, MAXFF) ; 
    if (nsplit ==0 ) continue ;
    sx = spt[0] ;
    if (sx[0] == '#')  { 
     freeup(spt, nsplit) ;
     continue ;
    }
    freeup(spt, nsplit) ;
    fclose(fff) ;
    if (nsplit>1) return NO ; 
    return YES ;
   }
   fatalx("(iseigenstrat) no genotyped data found\n") ;

}
int
ineigenstrat(char *gname, SNP **snpm, Indiv **indiv, int numsnps, int numind) 
// supports enhanced format fist character X => all missing data for SNP
{
   FILE *fff ;
   char *line, c ;
   char *spt[2], *sx ;
   int nsplit, rownum=0, k, num ;
   int maxstr, maxff = 2 ; 
   int nind, nsnp, len ;
   double y  ;
   unsigned char *buff  ;
   char *packit, *pbuff ;
   int *gtypes, g, g1, g2 ;
   SNP *cupt ; 
   Indiv *indx ;
   int nbad=0 ;


   packmode = YES ;
   maxstr = numind+10 ;
   ZALLOC(line, maxstr, char) ;

   nind = numind ; 
   nsnp = numsnps ;
   ZALLOC(gtypes, nind, int) ;

   y = (double) (nind * 2) / (8 * (double) sizeof (char)) ;  
   rlen = nnint(ceil(y)) ;
   rlen = MAX(rlen, 48)  ;  
   ZALLOC(buff, rlen, unsigned char) ;

   packlen = rlen*nsnp ;
   if (packgenos==NULL) {
    ZALLOC(packgenos, packlen, char) ;
    clearepath(packgenos) ;
   }

   openit(gname, &fff, "r") ;

   rownum = 0 ;
   pbuff = packgenos  ; 
   while (fgets(line, maxstr, fff) != NULL)  {
    nsplit = splitup(line, spt, maxff) ; 
    if (nsplit ==0 ) continue ;

    sx = spt[0] ;
    if (sx[0] == '#')  { 
     freeup(spt, nsplit) ;
     continue ;
    }


    if (nsplit>1) fatalx("(ineigenstrat) more than 1 field\n") ;


    if (rownum>=numsnps) fatalx("(ineigenstrat) too many lines in file %d %d\n", rownum, numsnps) ;
    num = pedcols[rownum] ;
    cupt = snpm[num] ; 
    ++rownum ;
    if (cupt == NULL) continue ;

     if (cupt -> ngtypes == 0) { 
      if (packmode == NO) {
       ZALLOC(cupt -> gtypes, numind, int) ;
       ivclear(cupt -> gtypes, -1, numind) ;
      }
      else  {
       ZALLOC(cupt -> gtypes, 1, int) ;
       cupt -> pbuff = pbuff ;  
       pbuff += rlen ;
      }
      cupt -> ngtypes = numind ;
     }

    if (sx[0] == 'X') { 
     freeup(spt, nsplit) ;
     continue ;
    }

    len = strlen(sx) ;  
    if (len != nind) {
       printf("(ineigenstrat) bad line %d ::%s\n", rownum, line) ;
       fatalx("(ineigenstrat) mismatch line length %d %d\n", len, nind) ;
    }

    for (k=0; k<len; k++) { 
     sscanf(sx+k, "%c", &c) ;  
     g = -2 ; 
     if (c=='0') g = 0 ;
     if (c=='1') g = 1 ;
     if (c=='2') g = 2 ;
     if (c=='9') g = -1 ;

     if (g==-2) fatalx("(ineigenstrat) bad character %c\n", c) ;

     if (indiv[k] -> ignore) g = -1 ;
     if (checkxval(cupt, indiv[k], g) == NO) g = -1 ;

    indx = indiv[k] ;
    if (checkxval(cupt, indx, g) == NO) g = -1 ;
     g2 = g ;
     if (g2<0) continue ;
     g1 = getgtypes(cupt, k) ;
     if ( (g1>=0) && (g1 != g2)) ++nbad ;
     putgtypes(cupt, k, g2) ;
    }
    freeup(spt, nsplit) ;
   }
   if (rownum != numsnps) fatalx("(ineigenstrat) mismatch in numsnps %d and numlines %d\n", numsnps, rownum) ;
   fclose(fff) ;
   return nbad ;
}

int calcishash(SNP **snpm, Indiv **indiv, int numsnps, int numind, int *pihash, int *pshash)
{
  char **arrx ;  
  int ihash, shash, n, num ;  
  int i ;
  Indiv *indx ; 
  SNP *cupt ;

   n = numind ;
   ZALLOC(arrx, n, char *) ;

   num = 0 ;
   for (i=0; i<n ; i++)  {  
    indx = indiv[i] ;  
    arrx[num] = strdup(indx -> ID) ;
    ++num ;
   }
   *pihash = hasharr(arrx, num) ;

   freeup(arrx, num) ;  
   free(arrx) ;
   
   n = numsnps ;
   ZALLOC(arrx, n, char *) ;
   num = 0 ;
   for (i=0; i<n ; i++)  {  
    cupt = snpm[i] ;  
    if (cupt -> isfake) continue ;
    arrx[num] = strdup(cupt -> ID) ;
    cupt -> ngtypes = numind ;
    ++num ;
   }
   *pshash = hasharr(arrx, num) ;
   freeup(arrx, num) ;  
   free(arrx) ;
   return num ;

}

void
inpack(char *gname, SNP **snpm, Indiv **indiv, int numsnps, int numind) 
{
   char **arrx, junk[10] ;  
   int n, num, ihash, shash, i, g, j, k, t ;
   int xihash, xshash, xnsnp, xnind ;
   int nind , nsnp, irec ;
   Indiv *indx ; 
   SNP *cupt ;
   double y ;
   unsigned char *buff  ;
   int fdes, ret ;
   char *packit, *pbuff ;

   nind = n = numind ;
   nsnp = calcishash(snpm, indiv, numsnps, numind, &ihash,  &shash) ; 

   y = (double) (nind * 2) / (8 * (double) sizeof (char)) ;  
   rlen = nnint(ceil(y)) ;
   rlen = MAX(rlen, 48)  ;  
   ZALLOC(buff, rlen, unsigned char) ;

  fdes = open(gname, O_RDONLY) ;
  if (fdes<0) {  
   perror("open failure") ;
   fatalx("(ispack) bad open %s\n", gname) ;
  } 
  t = read(fdes, buff, rlen ) ;
  if (t<0) {  
   perror("read failure") ;
   fatalx("(inpack) bad read") ;
  } 

  if (hashcheck) {
   sscanf((char *) buff,"GENO %d %d %x %x", &xnind, &xnsnp, &xihash, &xshash) ;
   if (xnind != nind) fatalx("OOPS number of individuals %d != %d in input files\n", nind, xnind) ;
   if (xnsnp != nsnp) fatalx("OOPS number of SNPs %d != %d in input file: %s\n", nsnp, xnsnp, gname) ;
   if (xihash != ihash) fatalx("OOPS indiv file has changed since genotype file was created\n") ;
   if (xshash != shash) fatalx("OOPS snp file has changed since genotype file was created\n") ;
  }

  packlen = rlen*nsnp ;
  ZALLOC(packgenos, packlen, char) ;
  clearepath(packgenos) ;

/**
   printf("packgenos: %x end: %x  len:  %d\n", 
    packgenos, packgenos+packlen-1, packlen) ;
*/

  t = read(fdes, packgenos, packlen ) ;
  if (t<0) {  
   perror("read failure") ;
   fatalx("(inpack) bad data read") ;
  } 
  if (t != packlen) {  
   perror("read failure (length mismatch)") ;
   fatalx("(inpack) bad data read (length mismatch) %d %d\n", t, packlen) ;
  } 
/* now set up pointers into packed data */
   pbuff = packgenos ;
   for (i=0; i<numsnps ; i++)  {  
    cupt = snpm[i] ;  
    if (cupt -> isfake) continue ;
    cupt -> pbuff = pbuff ;
    pbuff += rlen ;
// now check xhets
    for (k=0; k<numind; ++k) {  
     indx = indiv[k] ;  
     g = getgtypes(cupt, k) ;
     if (checkxval(cupt, indx, g) == NO) {
       putgtypes(cupt, k, -1) ;
     }
    }
   }

  free(buff) ;  
  close(fdes) ;
}

void
getsnpsc(char *snpscname, SNP **snpm, int numsnps)
{

   FILE *fff ;
   int  score  ;
   SNP *cupt ;
   char line[MAXSTR] ;
   char *spt[MAXFF], *sx ;
   int nsplit, num=0, k ;
   double y ;


   if (snpscname == NULL)  fatalx("no snpsc file\n") ;
   else openit(snpscname, &fff, "r") ;

   while (fgets(line, MAXSTR, fff) != NULL)  {
    nsplit = splitup(line, spt, MAXFF) ; 
    if (nsplit ==0 ) continue ;
    sx = spt[0] ;
    if (sx[0] == '#')  { 
     freeup(spt, nsplit) ;
     continue ;
    }
    k = snpindex(snpm, numsnps, sx) ;
    if (k<0) {  
     printf("*** warning.  snp %s in snpscname but not in main snp file\n", spt[0]) ;
     freeup(spt, nsplit) ;
     continue ;
    }
    y = atof(spt[1]) ;
    y += .1 * gauss() ; // dither 
    cupt = snpm[k] ; 
    cupt -> score = y ;
    freeup(spt, nsplit) ;
   }

   if (snpscname != NULL)  fclose(fff)  ;  
}

void setepath(SNP **snpm, int nsnps) 
{
 int i ; 
 SNP *cupt ;
 char *pbuff ;

 if (packlen==0) fatalx("(setepath) packlen unset\n") ;
 ZALLOC(packepath, packlen, char) ;
 printf("setepath. packlen: %d  rlen: %d\n", packlen, rlen) ;
 pbuff = packepath ;
 for (i=0 ; i<nsnps; i++)  {  
   cupt = snpm[i] ;  
   if (cupt -> isfake) continue ;
   cupt -> ebuff = pbuff ;
   pbuff += rlen ;
 }
 clearepath(packepath) ;
}
void clearepath(char *packp) 
{
 cclear((unsigned char *) packp, 0XFF, packlen) ;
}

int getpedgenos(char *gname, SNP **snpmarkers, Indiv **indivmarkers, 
 int numsnps, int numindivs, int nignore) 
{

 int val ;
 int ngenos = 0 ;

 SNP *cupt ;
 Indiv *indx ;

 char *line ;
 char **spt, *sx ;
 char  c ;
 int nsplit, num=0 ;
 int skipit ;
 int numf, snpnumber, nsnp ;
 int k, n, t, i ;
 FILE *fff ;
 int **gcounts, *gvar, *gref ;
 int xvar, xref ;
 int parity, colbase, ncols ; 
 int snpnum ;
 int markernum = -99 ;
 int n1, n2 ;

/**
 markernum = snpindex(snpmarkers, numsnps, "rs3002685") ;
 if (markernum <0) fatalx("qq1") ;
*/

  maxgenolinelength = MAX(maxgenolinelength, maxlinelength(gname)) ;

// printf("maxlinelen %d\n", maxlinelength(gname)) ;
 ZALLOC(line, maxgenolinelength+1, char) ;

 cleargdata(snpmarkers, numsnps, numindivs) ;
 nsnp = numsnps  ;

 ZALLOC(gcounts, nsnp, int *) ;
 for (i=0; i<nsnp; i++)  { 
   ZALLOC(gcounts[i], 5, int) ;
 }
 genopedcnt(gname, gcounts, nsnp) ;

 ZALLOC(gvar, nsnp, int) ;
 ZALLOC(gref, nsnp, int) ;
 setgref(gcounts, nsnp, gvar, gref) ;

 for (i=0; i<nsnp; ++i) { 
  cupt = snpmarkers[i] ;  
  if (cupt -> alleles[0] != CNULL) { 
   c =  cupt -> alleles[0] ;  
   gvar[i] = xpedval(c) ;
   c =  cupt -> alleles[1] ;  
   gref[i] = xpedval(c) ;
  }
  else {  
   c = x2base(gvar[i]) ;
   cupt -> alleles[0]  = c ;  
   c = x2base(gref[i]) ;
   cupt -> alleles[1]  = c ;  
  }
 }

 numf = 2*nsnp+10 ;
 ZALLOC(spt, numf, char *) ;

 openit(gname, &fff, "r") ;
 while (fgets(line, maxgenolinelength, fff) != NULL)  {
   nsplit = splitup(line, spt, numf) ; 
   if (nsplit == 0) continue ;
   sx = spt[0] ;
   skipit = NO ;
   if (sx[0] == '#') {  	
    freeup(spt, nsplit) ;
    continue ;
   }
   if (num==0) {  
    parity = nsplit % 2 ;  
    ncols = nsplit ; 
    colbase = 6 + parity ;
   }
   if (nsplit != ncols) 
    fatalx("bad number of fields %d %d\n", ncols, nsplit) ;

   for (k=colbase ; k < nsplit-1 ; k+=2)  { 
    snpnumber = (k-colbase)/2 ;

    if (snpnumber >= numpedcols) fatalx("pedcols overflow\n") ;
    snpnum = pedcols[snpnumber] ;
    if (snpnum<0) fatalx("logic bug (bad pedcols)\n") ;

    xvar = gvar[snpnum] ; 
    xref = gref[snpnum] ; 

    t = 0 ; 

    n1 = n = pedval(spt[k]) ; 
    n2 = pedval(spt[k+1]) ; 

    if ((n1==5) && (n2==5)) { 
     val = -1 ;
     putgtypes(cupt, num, val) ;
     continue ;
    }

    if ((n<0) || (n>4)) fatalx("(getpedgenos) %s bad geno %s\n", gname, spt[k]) ;
    if (n==xvar) ++t ; 
    if ((n != xvar) && (n != xref)) t = -10 ;

    n = n2 ;               
    if ((n<0) || (n>4)) fatalx("(getpedgenos) %s bad geno %s\n", gname, spt[k+1]) ;
    if (n==xvar) ++t ; 
    if ((n != xvar) && (n != xref)) t = -10 ;

    if (t<0) t = -1 ;
    cupt = snpmarkers[snpnum] ;  
    if (cupt -> ignore) continue ;
    val = t ;
    if (checkxval(cupt, indivmarkers[num], val) == NO) val = -1 ;
    putgtypes(cupt, num, val) ;
    if (val>=0) ++ngenos ;
   }
   freeup(spt, nsplit) ;
   ++num ;
  }

 free(spt) ;
 fclose(fff) ;

 for (i==0; i<nsnp; i++)  { 
   free(gcounts[i]) ;
 }

 free(gcounts) ;
 free(gref) ; 
 free(gvar) ;
 free(line) ;

 printf("genotype file processed\n") ;
 return ngenos ;
}

void
genopedcnt(char *gname, int **gcounts, int nsnp) 
{

 char *line ;
 char **spt, *sx ;
 int nsplit, num=0 ;
 int skipit ;
 int numf, snpnumber, snpnum ; 
 int k, n ;
 FILE *fff ;
 int parity, ncols, colbase ;  

// gcounts already zeroed 

  maxgenolinelength = MAX(maxgenolinelength, maxlinelength(gname)) ;

//  printf("maxlinelen %d\n", maxlinelength(gname)) ;
 ZALLOC(line, maxgenolinelength+1, char) ;

 numf = 2*nsnp+10 ;
 ZALLOC(spt, numf, char *) ;

 openit(gname, &fff, "r") ;
 while (fgets(line, maxgenolinelength, fff) != NULL)  {
   nsplit = splitup(line, spt, numf) ; 
   if (nsplit == 0) continue ;
   skipit = NO ;
   sx = spt[0] ;
   if (sx[0] == '#') { 
     freeup(spt, nsplit) ;
     continue ;
   }
   if (num==0) {  
    parity = nsplit % 2 ;  
    ncols = nsplit ; 
    colbase = 6 + parity ;
   }
   for (k=colbase ; k < nsplit-1 ; k+=2)  { 
    snpnumber = (k-colbase)/2 ;
    if (snpnumber >= numpedcols) fatalx("pedcols overflow\n") ;
    snpnum = pedcols[snpnumber] ;
    if (snpnum<0) fatalx("logic bug (bad pedcols)\n") ;
    n = pedval(spt[k]) ; 
//  if ((n<0) || (n>4)) fatalx("(genopedcnt) %s bad geno %s\n", gname, spt[k]) ;
    if ((n<0) || (n>4)) continue ;
    if (n>0) {
     ++gcounts[snpnum][n] ;
     ++num ;
    }
    n = pedval(spt[k+1]) ; 
//  if ((n<0) || (n>4)) fatalx("(genopedcnt) %s bad geno %s\n", gname, spt[k+1]) ;
    if ((n<0) || (n>4)) continue ;
    if (n>0) {
     ++gcounts[snpnum][n] ;
     ++num ;
    }
   }
   freeup(spt, nsplit) ;
   continue ;
 }
 free(spt) ;
 free(line) ;
 fclose(fff) ;
 return  ;
}

void
outfiles(char *snpname, char *indname, char *gname, SNP **snpm, 
  Indiv **indiv, int numsnps, int numindx, int packem, int ogmode) 
/** 
 call at end of main program usually 
*/
{
   int sizelimit = 10000000 ;
   int numind ;

   numind = rmindivs(snpm, numsnps, indiv, numindx) ;

   switch (outputmode)  {  

    case EIGENSTRAT: 
     printf("eigenstrat output\n") ;
     outeigenstrat(snpname, indname, gname, snpm, indiv,  numsnps, numind)  ;
     return ;

    case PED: 
     printf("ped output\n") ;
     outped(snpname, indname, gname, snpm, indiv,  numsnps, numind, ogmode)  ;
     return ;

    case PACKEDPED: 
     printf("packedped output\n") ;
     outpackped(snpname, indname, gname, snpm, indiv,  numsnps, numind, ogmode)  ;
     return ;

   case PACKEDANCESTRYMAP:
     if (snpname != NULL) printsnps(snpname, snpm, numsnps, indiv, NO, NO) ;
     packem = YES ;
     printdata(gname, indname, snpm, indiv, numsnps, numind, packem) ;
     return ;


    case ANCESTRYMAP:  
    default:  
     if (snpname != NULL) printsnps(snpname, snpm, numsnps, indiv, NO, NO) ;
     packem = NO ;
     if (numsnps > (sizelimit/numind)) packem = YES ;
     printdata(gname, indname, snpm, indiv, numsnps, numind, packem) ;
     return ;

   }

}

void
outeigenstrat(char *snpname, char *indname, char *gname, SNP **snpm, Indiv **indiv, int numsnps, int numind) 
{
   FILE *fff, *ifile  ;
   int  g, i, k ;
   SNP *cupt ; 
   Indiv *indx ;
   char ss[MAXSTR] ;


   settersemode(YES) ;
   if (snpname != NULL) 
    printsnps(snpname, snpm, numsnps, indiv, NO, NO) ;
  if (indname != NULL) {
   openit(indname, &ifile, "w") ;
   for(i = 0; i< numind; i++) {
      indx = indiv[i];
      if (indx->ignore) continue ;            
      strcpy(ss, indx -> egroup) ;
      if (qtmode) {  
        sprintf(ss, "%9.3f", indx -> rawqval) ;
      }
      fprintf(ifile,"%20s %c %10s",indx->ID, indx->gender,ss);
      fprintf(ifile,"\n") ;
      continue ;
   }
   fclose(ifile) ;
  }

   if (gname == NULL) return ;

   openit(gname, &fff, "w") ; 
   for (k=0; k< numsnps;  k++) {  
    cupt  = snpm[k] ;
     if (outputall==NO)  {
      if (ignoresnp(cupt)) continue ;
      if (cupt -> isrfake) continue ;
     }
    for (i=0; i<numind; i++) {
     indx = indiv[i] ; 
     if (indx->ignore) continue ;            
     g = getgtypes(cupt, i) ;
     if (g<0) g=9 ;  
     fprintf(fff, "%1d", g);
    }
    fprintf(fff, "\n") ;
   }
   fclose(fff) ;
}
void
setgref(int **gcounts, int nsnp, int *gvar, int *gref) 
{
   int tt[5] ;  
   int i, kmax ;

   for (i=0; i<nsnp; i++) {  
    copyiarr(gcounts[i], tt, 5) ;  
    tt[0] = -9999 ;  
    ivlmaxmin(tt, 5, &kmax, NULL) ; 
    gvar[i] = kmax ;  
    if (tt[kmax] == 0) gvar[i] = 5 ;
    tt[kmax] = -9999 ; 
    ivlmaxmin(tt, 5, &kmax, NULL) ; 
    gref[i] = kmax ;  
    if (tt[kmax] == 0) gref[i] = 5 ;
   }
}
void
cleargdata(SNP **snpmarkers, int numsnps, int numindivs) 
// wipe out all genotype data 
{
   int i , k ;  
   SNP *cupt ;
   char *pbuff ;
   double y ;


   y = (double) (numindivs * 2) / (8 * (double) sizeof (char)) ;  
   rlen = nnint(ceil(y)) ;
   rlen = MAX(rlen, 48)  ;  
   packlen = rlen*numsnps ;

   if (packlen <= 0) fatalx("bad packlen\n") ;

   if ((packmode) && (packgenos == NULL)) {
    ZALLOC(packgenos, packlen, char) ;
    clearepath(packgenos) ;
   }

   pbuff = packgenos ;


   for (i=0; i<numsnps; i++) {  
    cupt = snpmarkers[i] ;
//  if (cupt -> ignore) continue ;
    if (cupt -> ngtypes == 0) { 
     if (packmode == NO) {
      ZALLOC(cupt -> gtypes, numindivs, int) ;
     }
     else  {
      ZALLOC(cupt -> gtypes, 1, int) ;
      cupt -> pbuff = pbuff ;  
      pbuff += rlen ;
     }
     cupt -> ngtypes = numindivs ;
     for (k=0; k<numindivs; ++k) {
      putgtypes(cupt, k, -1) ;
     }
    }
   }
}

void
setgenotypename(char **gname, char *iname) 
{
 if (ispedfile(iname) == NO) return ;
 if ((*gname != NULL) && strcmp(*gname, "NULL") ==0)  {  
   *gname = NULL ;  
   return ;
 }
 if (*gname != NULL) return ; 
 *gname = strdup(iname) ;
}
int maxlinelength(char *fname) 
// linelength including \n 
{

  int len, maxlen ; 
  int nl, t ; 
  FILE *fff ;

  maxlen = -1 ; 

  len = 0 ; 
  nl = (int) (unsigned char) '\n' ;

  openit(fname, &fff, "r") ;
  while ((t = fgetc(fff)) != EOF) {  
   ++len ;
   if (t==nl) {  
    maxlen = MAX(maxlen, len) ;
    len = 0 ;
   }
  }
  return maxlen ;
}
void settersemode(int mode) 
{
 tersemode = mode ;
}

void outindped(char *indname, Indiv **indiv, int numind, int ogmode)  
{  
   FILE *fff, *ifile  ;
   int  g, i, k ;
   Indiv *indx ;
   char c ;
   int pgender, astatus ;
   int dcode = 1 ;

   if (indname == NULL) return ; 

    openit(indname, &ifile, "w") ;
    for(i = 0; i< numind; i++) {
      indx = indiv[i];
      if (indx->ignore) continue ;            
      fprintf(ifile, "%6d %12s", i+1, indx->ID) ;
      fprintf(ifile, " %d %d", 0, 0) ;  // parents 
      c = indx->gender ;
      pgender = 0 ; 
      if (c == 'M') pgender = 1 ;
      if (c == 'F') pgender = 2 ;
      fprintf(ifile, " %d", pgender) ;
      if (ogmode == NO) {
       astatus = indx -> affstatus + 1 ;
       if (qtmode) {  
        fprintf(ifile, "%9.3f", indx -> rawqval) ;
       }
       else fprintf(ifile, " %d", astatus) ;
      }
      if (ogmode == YES) fprintf(ifile, " %10s", indx -> egroup) ;
      if (sevencolumnped) fprintf(ifile, " %d", dcode) ;
      fprintf(ifile, "\n") ;
    }
    fclose(ifile) ;
}

void
outped(char *snpname, char *indname, char *gname, SNP **snpm, Indiv **indiv, 
  int numsnps, int numind, int ogmode) 
{

   FILE *fff, *ifile  ;
   int  g, i, k ;
   SNP *cupt ; 
   Indiv *indx ;
   char c ;
   int pgender, astatus ;
   int g1, g2, dcode=1  ;

   settersemode(YES) ;
   if (snpname != NULL) 
    printmap(snpname, snpm, numsnps, indiv) ;
   if (indname!=NULL)
   outindped(indname, indiv, numind, ogmode) ;

   if (gname == NULL) return ;
    openit(gname, &fff, "w") ;
    for(i = 0; i< numind; i++) {
      indx = indiv[i];
      if (indx->ignore) continue ;            
      fprintf(fff, "%6d %12s", i+1, indx->ID) ;
      fprintf(fff, " %d %d", 0, 0) ;  // parents 
      c = indx->gender ;
      pgender = 0 ; 
      if (c == 'M') pgender = 1 ;
      if (c == 'F') pgender = 2 ;
      fprintf(fff, " %d", pgender) ;
      if (ogmode == NO) {
       astatus = indx -> affstatus +1 ;
       if (qtmode) {  
        fprintf(fff, "%9.3f", indx -> rawqval) ;
       }
       else fprintf(fff, " %d", astatus) ;
      }
      if (ogmode == YES) fprintf(fff, " %10s", indx -> egroup) ;
      if (sevencolumnped) fprintf(fff, " %d", dcode) ;
      for (k=0; k<numsnps; k++) { 
       cupt = snpm[k] ; 
      if (outputall==NO)  {
       if (ignoresnp(cupt)) continue ;
       if (cupt -> isrfake) continue ;
      }
       g = getgtypes(cupt, i) ;
       gtox(g, cupt -> alleles, &g1, &g2) ;
       fprintf(fff, "  %d %d", g1, g2 ) ;
       if ((g1>4) || (g2>4))  {  
        fprintf(fff, "\n") ;
        fflush(fff) ;
        fclose(fff) ;
        printf("bad genotype for snp %s alleles: ", cupt -> ID) ;
        printalleles(cupt, stdout) ;
        printf("\n") ;
        fatalx("trying to make invalid ped file %s\n", gname) ;
       }
      }
      fprintf(fff, "\n") ;
    }
    fclose(fff) ;
}

void gtox(int g, char *cvals, int *p1, int *p2) 
// output values for ped file using allele array
{

       int g1, g2 ;

       
       switch (g)  {
        case -1:
         *p1 = *p2 = 0; return ;
        case 0:
         g1 = 1 ; g2 = 1; break ;
        case 1:
         g1 = 1 ; g2 = 2; break ;
        case 2:
         g1 = 2 ; g2 = 2; break ;
        default:
         fatalx("(outped) bug %d\n", g) ;
       }
       if (cvals != NULL)  {  
        g1 = 3-g1 ; 
        g2 = 3-g2 ;
        g1 =  xpedval(cvals[g1-1]) ;
        g2 =  xpedval(cvals[g2-1]) ;
       }

       *p1 = MIN(g1, g2) ;
       *p2 = MAX(g1, g2) ;

}

void
outpackped(char *snpname, char *indname, char *gname, SNP **snpm, Indiv **indiv, 
  int numsnps, int numind, int ogmode) 
{

   FILE *fff, *ifile  ;
   int  g, i, k ;
   SNP *cupt ; 
   Indiv *indx ;
   char c ;
   int pgender, astatus ;
   int g1, g2, dcode=1  ;
   unsigned char ibuff[3]  ;
   unsigned char *buff ;
   int fdes, ret, blen ;
   int *gtypes ;
   double y ;


   settersemode(YES) ;
   if (snpname!=NULL)
    printmap(snpname, snpm, numsnps, indiv) ;
   if (indname != NULL)  
   outindped(indname, indiv, numind, ogmode) ;

   if (gname == NULL) return ;

   ibuff[0] = 0x6C ; 
   ibuff[1] = 0x1B ;
   ibuff[2] = 0x01 ;  

/** magic constants for snp major bed file */

   y = (double) (numind * 2) / (8 * (double) sizeof (char)) ;  
   blen = nnint(ceil(y)) ;
   ZALLOC(buff, blen, unsigned char) ;
   ZALLOC(gtypes, numind, int) ;

   fdes = open(gname, O_CREAT | O_TRUNC | O_RDWR, 0666);
 
   if (fdes<0) {
    perror("bad gname") ;
    fatalx("open failed for %s\n", gname) ;
   }
   if (verbose) 
    printf("file %s opened\n", gname) ;
 
   if (fdes<0) {
    perror("bad genoout") ;
    fatalx("open failed for %s\n", gname) ;
   }

   if (verbose) 
    printf("file %s opened\n", gname) ;

   ret = write(fdes, ibuff, 3) ;

   for (i=0; i<numsnps; i++) {  
    cupt = snpm[i] ;
      if (outputall==NO)  {
       if (ignoresnp(cupt)) continue ;
       if (cupt -> isrfake) continue ;
      }
    for (k=0; k<numind; ++k) { 
     g = getgtypes(cupt, k) ;
     if (g>=0) g = 2 - g ; 
     gtypes[k] = g ;
    }
    setbedbuff((char *) buff, gtypes, numind) ;
    ret = write(fdes, buff, blen) ;
     if (ret<0) {  
      perror("write failure") ;
      fatalx("(outpackped) bad write") ;
     } 
   }


   free(buff) ;
   close(fdes) ;
}
void
setbedbuff(char *buff, int *gtypes, int numind )
{
   int i, k ;  
   double y ;  
   int blen, wnum, wplace, bplace, t, g ;
   unsigned char c ;

   y = (double) (numind * 2) / (8 * (double) sizeof (char)) ;  
   blen = nnint(ceil(y)) ;

   c = 0xAA  ;  // missing 
   cclear((unsigned char *) buff, c, blen) ;

   for (k=0; k<numind; k++) {  
    wnum = k/4 ;
    t = k%4 ;
    wplace = 3-t ;  // switch for bed 
    bplace = 4*wnum + wplace ;
    g = bedval(gtypes[k]) ;
    wbuff((unsigned char *) buff, bplace, g) ;
   }
}
int bedval(int g) 
{
 if (g<0) return  1  ;  
 if (g==2) return 3  ;
 if (g==1) return 2 ;
 if (g==0) return 0 ;

 fatalx("(bedval) bad g value %d\n", g) ;

}

void
atopchrom(char *ss, int chrom) 
// ancestry chromosome -> map covention  
{
  switch (chrom) {  

   case 23: 
    strcpy(ss, "X") ;
    return ;

   case 24: 
    strcpy(ss, "Y") ;
    return ;

   default: 
    sprintf(ss, "%d", chrom) ;
  }
}
int
ptoachrom(char *ss) 
// map -> ancestry  
{
  char c ;

  c = ss[0] ; 

  if (c=='X') return 23 ;
  if (c=='Y') return 24 ;
  return atoi(ss) ;
}


void
printmap(char *snpname, SNP **snpm, int numsnps, Indiv **indiv) 
{ 

  char ss[5] ;
  int i ;  
  FILE *fff ;
  SNP *cupt ;
  char  c ;

  if (snpname == NULL) return  ;
  openit(snpname, &fff, "w") ;
  for (i=0; i<numsnps; i++) {  
   cupt = snpm[i] ;
      if (outputall==NO)  {
       if (ignoresnp(cupt)) continue ;
       if (cupt -> isrfake) continue ;
      }
   atopchrom(ss, cupt -> chrom) ;
   fprintf(fff, "%-2s", ss) ;
   fprintf(fff, " %12s", cupt -> ID) ;
   fprintf(fff, " %12.6f", cupt -> genpos) ;
   fprintf(fff, " %12.0f", cupt -> physpos) ;
   printalleles(cupt, fff) ;
   fprintf(fff, "\n") ;
  }
  fclose(fff) ;
}
char x2base(int x) 
// 12345 -> ACGTX
{
 char *blist = "?ACGT" ;
 if (x<0) return '?' ;
 if (x>4) return 'X' ;
  return blist[x] ;
}

int xpedval(char c) 
{
  char bb[2] ; 

  bb[1] = '\0' ;
  bb[0] = c ;

  if (isdigit(c)) return atoi(bb) ;
  return pedval(bb) ;



}
int pedval(char *sx) 
{
  char c ;  

  c = sx[0] ;
  if (c=='A') return 1 ;
  if (c=='C') return 2 ;
  if (c=='G') return 3 ;
  if (c=='T') return 4 ;
  if (c=='X') return 5 ;
  if (c=='N') return 5 ;

  if (c=='1') return 1 ;
  if (c=='2') return 2 ;
  if (c=='3') return 3 ;
  if (c=='4') return 4 ;
  if (c=='0') return 0 ;

  return 9  ;

}

int getbedgenos(char *gname, SNP **snpmarkers, Indiv **indivmarkers, 
 int numsnps, int numindivs, int nignore) 
{

 int val, i, k, x ;
 int t, wnum, wplace ;
 int nsnp  ;
 int ngenos = 0 ;

 SNP *cupt ;
 Indiv *indx ;

 unsigned char *buff, ibuff[3], jbuff[3] ;
 double y ;  
 int blen ;
 int fdes ;


  ibuff[0] = 0x6C ; 
  ibuff[1] = 0x1B ;
  ibuff[2] = 0x01 ;  


 cleargdata(snpmarkers, numsnps, numindivs) ;
 nsnp = numsnps  ;

  y = (double) (numindivs * 2) / (8 * (double) sizeof (char)) ;  
  blen = nnint(ceil(y)) ;

  ZALLOC(buff, blen, unsigned char) ;

  fdes = open(gname, O_RDONLY) ;

  if (fdes<0) {  
   perror("open failure") ;
   fatalx("(getbedgenos) bad open %s\n", gname) ;
  } 

  t = read(fdes, jbuff, 3 ) ;

  if (t<0) {  
   perror("read failure") ;
   fatalx("(getbedgenos) bad read") ;
  } 


// check magic
 for (k=0; k<3; k++) {
  if (ibuff[k] != jbuff[k]) {
   fprintf(stderr, "magic failure: ") ;
   fprintf(stderr, " %x %x %x", jbuff[0], jbuff[1], jbuff[2]) ;
   fprintf(stderr, " %x %x %x", ibuff[0], ibuff[1], ibuff[2]) ;
   fprintf(stderr, "\n") ;
   fatalx("(getbedgenos) magic failure\n") ;
  }
 }



  for (i=0; i<nsnp; i++) {  
   cupt = snpmarkers[i] ;
   t = read(fdes, buff, blen) ; 

   if (t<0) {  
    perror("read failure") ;
    fatalx("(getbedgenos) bad read") ;
   } 
   if (cupt -> ignore) continue ;

   for (k=0; k<numindivs; k++) {  
    indx = indivmarkers[k] ;
    wnum = k/4 ;
    t = k%4 ;
    wplace = 3-t ;  // switch for bed 
    wplace += 4*wnum ;
    x = rbuff(buff, wplace) ;
    val = ancval(x) ;
    if (checkxval(cupt, indx, val) == NO) val = -1 ;
    putgtypes(cupt, k, val) ;
    if (val >= 0) ++ngenos ;
   }
  }

 free(buff) ;
 printf("genotype file processed\n") ;
 return ngenos ;
}

int ancval(int x) 
// bed -> anc 
// 1/22/07 allele flipped
{
 if (x==1) return -1 ;
 if (x==3) return 0 ;
 if (x==2) return 1 ;
 if (x==0) return 2 ;
 fatalx("(ancval) bad value %d\n", x) ;
}



void setomode(enum outputmodetype *outmode, char *omode) 
{
 char *ss ;
 int len, i ;
 if (outmode == NULL) return ; 
 *outmode = PACKEDANCESTRYMAP ;
 if (omode == NULL) return ;

 ss = strdup(omode) ;
 len = strlen(ss) ;
 for (i=0; i<len ; i++) {
  ss[i] = tolower(ss[i]) ;
 }

 if (strcmp(ss, "eigenstrat") == 0)  *outmode = EIGENSTRAT ;
 if (strcmp(ss, "ascii") == 0)  *outmode = EIGENSTRAT ;
 if (strcmp(ss, "alkes") == 0)  *outmode = EIGENSTRAT ;
 if (strcmp(ss, "ped") == 0)  *outmode = PED ;
 if (strcmp(ss, "packedped") == 0)  *outmode = PACKEDPED ;
 if (strcmp(ss, "packedancestrymap") == 0)  *outmode = PACKEDANCESTRYMAP ;
 if (strcmp(ss, "ancestrymap") == 0)  *outmode = ANCESTRYMAP ;

 free(ss) ;
}
void snpdecimate(SNP **snpm, int nsnp, int decim, int mindis, int maxdis)  
{
 int chrom = -1 ; 
 SNP **cbuff, *cupt, *cupt2 ; 
 int k, k2, n, t  ;
  
 printf( "snpdecimate called: decim: %d  mindis: %d  maxdis: %d\n", decim, mindis, maxdis) ; 
 ZALLOC(cbuff, nsnp, SNP *) ;
 for (k=0; k<nsnp; ++k)  {  
  cupt = snpm[k] ;    
  if (cupt -> chrom != chrom)  {  
   chrom = cupt -> chrom ;
   n = 0 ;
   for (k2=k; k2 <nsnp; ++k2)   {  
    cupt2 = snpm[k2] ;
    if (cupt2 -> chrom != chrom) break ;
    if (cupt2 -> ignore) continue ;  
    if (cupt2 -> isfake) continue ;
    cbuff[n] = cupt2 ;  
    ++n ;
   }
   if (n<decim) continue ;
   decimate(cbuff, n, decim, mindis, maxdis) ;
  }
 }
}
void 
decimate(SNP **cbuff, int n, int decim, int mindis, int maxdis) 
{ 
 int k, t, u, dis, len ; 
 int *ttt ;
 SNP *cupt ; 

 cupt = cbuff[0] ;
 if (n<2) return ;
 if (n<decim) return  ;
 ZALLOC(ttt, n, int) ;
 for (k=1; k<n; ++k)  {  
  dis = (int) (cbuff[k] -> physpos - cbuff[k-1] -> physpos) ; 
  if (dis > maxdis)  {  
   decimate(cbuff, k-1, decim, mindis, maxdis) ;
   decimate(cbuff+k, n-k, decim, mindis, maxdis) ;
   return ;
  }
 }
 t = ranmod(decim) ;
 ttt[t] = 1 ;

 u = t + decim ; 

 for (;;) {  
  if (u>=n) break ;
  dis = (int) (cbuff[u] -> physpos - cbuff[t] -> physpos) ; 
   if (dis<mindis) { 
    ++u ;
    continue ; 
   }
   len = u-t-1 ;
   ivclear(ttt+t+1, 1, len) ;  
   t = u ;
   u = t + decim ;
 }
 for (k=0 ; k < n; ++k)  {  
  if (ttt[k] == 1) cbuff[k] -> ignore = YES ;
 }
// debug 
 if (verbose) {
  for (k=0 ; k < n; ++k)  {  
   printf("zz %6d %20s %20d  %3d\n", k, cbuff[k] -> ID, (int) cbuff[k] -> physpos, ttt[k]) ;
  }
 }

 free(ttt) ;

}

int killhir2(SNP **snpm, int numsnps, int numind, double physlim, double genlim, double rhothresh)  
/** 
 physlim = genlim = 0 => kill monomorphs 
*/
{
 double *badbuff ; 
 int  *xbadbuff ;
 SNP *cupt, *cupt1, *cupt2 ;

#define BADBUFFSIZE 100000 ;
 int badbuffsize = BADBUFFSIZE ;
 int i,j, k, nbad, kmax, kmin, t, j1, j2, lo, hi  ;
 int *gtypes ;
 double *x1, *x2, mean, dis, *p1 ;
 int nkill = 0, tj ;
 double y1, y2, y, rho, smax ;
 double **xx1, *yy1 ;
 SNP **snpxl ;

 if (physlim<0) return 0 ;
 if (genlim<0) return 0 ;

/** 
 step 1 give score to each SNP 
*/
 for (i=0; i<numsnps; i++) {  
   cupt = snpm[i] ; 
   if (cupt -> ignore) continue ; 
   cupt -> score = numvalidgtypes(cupt)  ;  
   cupt -> score += DRAND() ;  // jitter
 }

 ZALLOC(badbuff, badbuffsize, double) ; 
 ZALLOC(xbadbuff, badbuffsize, int) ; 
 ZALLOC(x1, numind, double) ;
 ZALLOC(x2, numind, double) ;
 ZALLOC(gtypes, numind, int) ;

 xx1 = initarray_2Ddouble(10000, numind, 0.0) ;
 ZALLOC(yy1, 10000, double) ;
 ZALLOC(snpxl, 10000, SNP *) ;

 for (i=0; i<numsnps; i+=5000)  { 


  lo = i ;  
  hi = i+10000-1 ;
  hi = MIN(hi, numsnps-1) ;

  for (j=lo; j <=hi ; ++j)   {  
    p1 = xx1[j-lo] ;
    cupt = snpm[j] ; 
    snpxl[j-lo] = cupt ;
    grabgtypes(gtypes,  cupt,  numind) ; 
    floatit(p1, gtypes, numind) ;
    vvadjust(p1, numind, NULL) ; 
    y1 = asum2(p1, numind) ;  
    yy1[j-lo] = y1 ;
    if (y1<0.01) {  
     ++nkill ;
     cupt -> ignore = YES ;
    }
  }
  for (j1 = 0 ; j1 < 5000; ++j1) {  
   if (j1>(hi-lo)) break ;
   cupt1 = snpxl[j1] ;
   if (cupt1 -> ignore) continue ;
   nbad = 0 ;
   tj = 0 ;
   for (j2=j1+1; j2 <= hi-lo; ++j2) {  
    cupt2 = snpxl[j2] ;
    if (cupt2 -> ignore) continue ;
    if (cupt2 -> chrom != cupt1 -> chrom) break ;

    dis = cupt2 -> genpos - cupt1 -> genpos ;
    if (dis > genlim) break ; 

    dis = cupt2 -> physpos - cupt1 -> physpos ;
    if (dis > physlim) break ; 
    ++tj ;

    y1 = yy1[j1] ;
    y2 = yy1[j2] ;
    

    y = vdot(xx1[j1], xx1[j2], numind) / sqrt(y1*y2) ;  
    rho = y * y ;
    if (rho < rhothresh) continue ;
    badbuff[nbad] = cupt2 -> score ;  
    xbadbuff[nbad] = j2+lo ;  
    ++nbad ;

   }
   t = (j1+lo) % 100 ;  
   if (nbad == 0) continue ;  
   vlmaxmin(badbuff, nbad, &kmax, &kmin) ;  
   smax = snpm[kmax] -> score ; 
   if (smax > cupt1 -> score)  {  
    cupt1 -> ignore = YES ;  
    ++nkill ;
    continue ;
   }
   for (k=0; k<nbad; ++k)  { 
    j = xbadbuff[k] ;
    snpm[j] -> ignore = YES ;
    ++nkill ;
   }
  }
 }


 free2D(&xx1, 10000) ; 
 free(yy1) ; 
 free(snpxl) ;
 free(gtypes) ;
 free(badbuff) ; 
 free(xbadbuff) ;
 free(x1) ;
 free(x2) ;

 for (i=0; i<numsnps; i++) { 
  cupt = snpm[i] ;
  cupt -> score = 0.0 ;
 }

 printf("killr2 complete\n") ;
 return nkill ;
}

int vvadjust(double *cc, int n, double *pmean) 
/* take off mean  force missing to zero */
/* simpler version of vadjust */
{
 double ynum, ysum, y, ymean ;  
 int i, nmiss=0 ;

 ynum = ysum = 0.0 ;
 for (i=0; i<n; i++) {  
  y = cc[i] ;
  if (y < 0.0) { 
   ++nmiss ;
   continue ;
  }
  ++ynum ;
  ysum += y ;
 }
 if (ynum <= 1.5) {  
/* no data or monomorphic */
  vzero(cc, n) ; 
  if (pmean != NULL) *pmean = ysum/(ynum+1.0e-8) ;
  return nmiss ;
 }
 ymean = ysum/ynum ;
 for (i=0; i<n; i++) {  
  y = cc[i] ;
  if (y < 0.0) cc[i] = 0.0 ;
  else cc[i] -= ymean ;
 }
 if (pmean != NULL) *pmean = ymean ;
 return nmiss ;
}
static int setskipit(char *sx) 
{
   int skipit = NO ;
   if (sx[0] == '#') skipit = YES ;
   if (strcmp(sx,"SNP_ID") == 0) skipit = YES ;
   if (strcmp(sx,"Indiv_ID") == 0) skipit = YES ;
   if (strcmp(sx,"Chr") == 0) skipit = YES ;
   return skipit ;
}

int
inpack2(char *gname, SNP **snpm, Indiv **indiv, int numsnps, int numind) 
// load up packed genotype file for merge.
{
   char **arrx, junk[10] ;  
   int n, num, ihash, shash, i, g, j, k, t, g1, g2 ;
   int xihash, xshash, xnsnp, xnind ;
   int nind , nsnp, irec ;
   Indiv *indx ; 
   SNP *cupt, *cupt2 ;
   SNP xsnp ;
   double y ;
   unsigned char *buff, *tbuff  ;
   int fdes, ret ;
   char *packit, *pbuff ;
   int nbad = 0 ;
   n = numind ;

   ZALLOC(arrx, n, char *) ;

   num = 0 ;
   for (i=0; i<n ; i++)  {  
    indx = indiv[i] ;  
    arrx[num] = strdup(indx -> ID) ;
    ++num ;
   }
   ihash = hasharr(arrx, num) ;
   nind= num ;


   freeup(arrx, num) ;  
   free(arrx) ;
   
   n = numsnps ;
   ZALLOC(arrx, n, char *) ;
   num = 0 ;
   for (i=0; i<n ; i++)  {  
    cupt = snpm[i] ;
    arrx[num] = strdup(cupt -> ID) ;
    ++num ;
   }
   shash = hasharr(arrx, num) ;
   nsnp = num ;
   freeup(arrx, num) ;  
   free(arrx) ;

   y = (double) (nind * 2) / (8 * (double) sizeof (char)) ;  
   rlen = nnint(ceil(y)) ;
   rlen = MAX(rlen, 48)  ;  
   ZALLOC(buff, rlen, unsigned char) ;
   ZALLOC(tbuff, rlen, unsigned char) ;

  fdes = open(gname, O_RDONLY) ;
  if (fdes<0) {  
   perror("open failure") ;
   fatalx("(inpack2) bad open %s\n", gname) ;
  } 
  t = read(fdes, buff, rlen ) ;
  if (t<0) {  
   perror("read failure") ;
   fatalx("(inpack2) bad read") ;
  } 
  sscanf((char *) buff,"GENO %d %d %x %x", &xnind, &xnsnp, &xihash, &xshash) ;

  if (xnind != nind) fatalx("(inpack2) nind mismatch %d %d \n", nind, xnind) ;
  if (xnsnp != nsnp) fatalx("(inpack2) nsnp mismatch\n") ;
  if (xihash != ihash) fatalx("(inpack2) ihash mismatch\n") ;
  if (xshash != shash) fatalx("(inpack2) shash mismatch\n") ;


/* now copy genotypes */
   for (i=0; i<n ; i++)  {  
    t = read(fdes, tbuff, rlen) ;
    if (t != rlen) {  
     perror("read failure") ;
     fatalx("(inpack2) bad data read") ;
    }  
    cupt = snpm[i] ;  
    if (cupt -> isfake) continue ;
    xsnp = *cupt ;
    cupt2 = &xsnp ;
    cupt2 -> pbuff = (char *) tbuff ;
    for (k=0; k<numind; ++k) {  
     g2 = getgtypes(cupt2, k) ;
     if (g2<0) continue ;
     g1 = getgtypes(cupt, k) ;
     if ( (g1>=0) && (g1 != g2)) ++nbad ;
     putgtypes(cupt, k, g2) ;
    }

// now check xhets
    for (k=0; k<numind; ++k) {  
     if (cupt -> chrom != 23) break; 
     indx = indiv[k] ;  
     g = getgtypes(cupt, k) ;
     if (checkxval(cupt, indx, g) == NO) {
       putgtypes(cupt, k, -1) ;
     }
    }
   }

  free(buff) ;  
  free(tbuff) ;  
  close(fdes) ;
  return nbad ;
}

void
getgenos_list(char *genotypelist, SNP **snpmarkers, Indiv **indivmarkers, 
 int  numsnps, int numindivs, int nignore) 
{
   char **fnames, *fn ;
   int n ; 
   int k, nbad, isok ; 

   dofreeped = NO ;
   n = numlines(genotypelist) ;
   ZALLOC(fnames, n, char *) ;
   n = getlist(genotypelist, fnames) ;
   getgenos(fnames[0], snpmarkers, indivmarkers, 
     numsnps, numindivs, nignore) ;

   for (k=1; k<n; ++k) {  
    fn = fnames[k] ;
    isok = NO ;
    if (ispack(fn))  {
     nbad = inpack2(fn, snpmarkers, indivmarkers, numsnps, numindivs)  ;          
     isok = YES ;
    }
    if (iseigenstrat(fn))  {
     nbad = ineigenstrat(fn, snpmarkers, indivmarkers,  numsnps, numindivs)  ; 
     isok = YES ;
    }
    if (nbad>0) printf("%s genotypes mismatches: %d\n", fn, nbad) ;
    if (isok == NO) fatalx("file %s must be packed or eigenstrat format\n") ;
   }

   dofreeped = YES ;
   freeped() ;

   free(fnames) ;
}

int str2chrom(char *ss) 
{
  if (strcmp(ss, "X") == 0) return 23 ;
  if (strcmp(ss, "Y") == 0) return 24 ;
  if (strcmp(ss, "MT") == 0) return MTCHROM ;
  if (strcmp(ss, "XY") == 0) return XYCHROM ;
  if (!isnumword(ss)) return -1 ;
  return atoi(ss) ;
}

int checksize(int numindivs, int numsnps, enum outputmodetype outputmode) 
{
// -1 try packed format

  double y ;
  long z ;

  if (sizeof(z) == 8) checksizemode = NO ;
  if (checksizemode == NO) return 1 ; 
  
  y = (double) numindivs ; 
  y *= (double) numsnps  ;
  
  if (y>8.0e9) return -2 ;  

  switch (outputmode)   { 

   case ANCESTRYMAP:  
    if (y>5.0e7)  return -1 ;
    break ;
   case EIGENSTRAT:  
    if (y>2.0e9) return -1 ;
    break ;
   case PED:  
    if (y>4.0e8) return -1 ;
    break ;
   case PACKEDPED:  
    break ;
   case PACKEDANCESTRYMAP:  
    break ;
   default:
    fatalx("unknown outputmode\n") ;
  }
  return 1 ;
}

int snprawindex(SNPDATA **snpraw, int nreal, char *sname) 
{
     int k ; 
     char **ss ;  

     freesnpindex() ;

     if (snprawtab==NO) {  
// set up hash table
      snprawtab = YES ;
      ZALLOC(ss, nreal, char *) ;
      for (k=0; k< nreal; k++) {  
       ss[k] = strdup(snpraw[k] -> ID) ;
      }
      xloadsearch(ss, nreal) ;
      freeup(ss, nreal) ;
      free(ss) ;
     }
     
     k = xfindit(sname) ;
     return k ;
}
void freesnprawindex() 
{
 if (snprawtab == NO) return ;
 snprawtab = NO ; 
 xdestroy() ;
}

void cntpops(int *count, Indiv **indm, int numindivs, char **eglist, int numeg)
// count number of samples for each pop
{
    Indiv *indx ;
    int t, k ;

    ivzero(count, numeg) ;
    for (k=0; k<numindivs; ++k) { 
     indx = indm[k] ; 
     if (indx -> ignore) continue ;
     t = indxindex(eglist, numeg, indx -> egroup) ;  
     if (t<0) continue ;
     ++count[t] ;
    }
}

char *getpackgenos() 
{
 return packgenos ;
}

void clearpackgenos() 
{
 packgenos  = NULL ;
}


void genocloseit(genofile *gfile) 
{

 genofile *gpt ;
 SNP *cupt ; 
 int i ;
 gpt = gfile ;
 
 free(gpt -> buff) ;
 for (i=0; i< gpt -> numsnps; i++) { 
  cupt = gpt -> snpm[i] ;
  freecupt(&cupt) ;  
 }
 free(gpt -> snpm) ;

 for (i=0; i< gpt -> numindivs; i++) { 
  free(gpt -> indivm[i]) ;
 }
 free(gpt -> indivm) ;

 close(gpt -> fdes) ;

}

int 
genoopenit(genofile **gfile, char *geno2name, SNP **snp2m, Indiv **indiv2m, int numsnp2, int numindiv2, int nignore) 
/* only one gfile can be open */
{
 static genofile xfile  ;
 genofile *gpt ;
 double  y ;
 int rlen, fdes, t ;
 static unsigned char *buff ;
 int xihash, xshash, xnsnp, xnind ;
 int ihash, shash ;
 char *gname  ;
 int nsnp, nind ;
 

 if (geno2name == NULL) fatalx("(genoopenit) null name\n") ;
 if (!ispack(geno2name)) fatalx("(genoopenit) not packed ancestrymap format\n") ;
 gpt = *gfile = &xfile ;
 strcpy(gpt -> gname, geno2name) ; 
 gpt -> snpm = snp2m ;
 gpt -> indivm = indiv2m  ;
 gpt -> numsnps = numsnp2 ;
 gpt -> numindivs = numindiv2 ;
 

 y = (double) (numindiv2 * 2) / (8 * (double) sizeof (char)) ;
 rlen = nnint(ceil(y)) ;

 gpt -> rlen = rlen ;
 rlen = MAX(rlen, 48) ;
 ZALLOC(buff, rlen, unsigned char) ;
 gpt -> buff = buff ;

 fdes = open(geno2name, O_RDONLY) ;  
 if (fdes<0) return fdes ;
 gpt -> fdes = fdes ;
 gpt -> snpindex  = -1 ;

 t = read(fdes, buff, rlen) ;  
 if (t<0) fatalx("(genoopenit) bad initial read\n") ;

 nsnp = numsnp2 ;
 nind = numindiv2 ;
 gname = geno2name ;

 calcishash(snp2m, indiv2m, nsnp, nind, &ihash,  &shash) ; 
 if (hashcheck)  {
  sscanf((char *) buff,"GENO %d %d %x %x", &xnind, &xnsnp, &xihash, &xshash) ;
  if (xnind != nind) fatalx("OOPS number of individuals %d != %d in input files\n", nind, xnind) ;
  if (xnsnp != nsnp) fatalx("OOPS number of SNPs %d != %d in input file: %s\n", nsnp, xnsnp, gname) ;
  if (xihash != ihash) fatalx("OOPS indiv file has changed since genotype file was created\n") ;
  if (xshash != shash) fatalx("OOPS snp file has changed since genotype file was created\n") ;
 }

 return 0 ;
/** 

typedef struct {  
 char gname[IDSIZE] ;  
 SNP **snpm ;
 Indiv **indivm ;  
 int numsnps; 
 int numindivs ; 
 int rlen ;
 int fdes ; 
 unsigned char *buff ;
 int snpindex ;
} genofile ;
*/

}

int 
genoreadit(genofile *gfile, SNP **pcupt)  
/** 
 return code 
 < 0 bad read  
 0 EOF 
 rlen good read 
*/
{
 genofile *gpt ;
 SNP *cupt ;
 int t, rlen, snum ;
 int k ;
 
 cupt = *pcupt = NULL ;
 gpt = gfile ;
 rlen = gpt -> rlen ; 
 t = read(gpt -> fdes, gpt -> buff, rlen) ;  
 if (t<0) fatalx("(genoreadit) bad read \n") ;
 if (t==0) return  0  ; 
 if (t< gpt -> rlen) fatalx("(genoopenit) premature EOF\n") ;
 ++gpt -> snpindex ;
 snum = gpt -> snpindex ;
 cupt = *pcupt = gpt -> snpm[snum] ;
 cupt -> tagnumber = snum ;
 cupt -> pbuff = (char *) gpt -> buff ;
 cupt -> ngtypes = gpt -> numindivs ;
 if (cupt -> gtypes == NULL) ZALLOC(cupt -> gtypes, 1, int) ;
 return rlen ;
}             


void putped(int num) 
{
  int *pp ;  
  int t ;

  pp = pedcolsa[num] ; 
  if (pp != NULL) free(pp) ;
  pp = NULL ;
  t = numpedcolsa[num] = numpedcols ;
  if (t==0) return ; 
  ZALLOC(pedcolsa[num], t, int) ;
  pp = pedcolsa[num] ;
  copyiarr(pedcols, pp, t) ;
}
void getped(int num) 
{
  int *pp ;  
  int t ;
  
  pp = pedcols ;
  if (pp != NULL) free(pp) ;
  pp = NULL ;
  t = numpedcols = numpedcolsa[num]  ;
  if (t==0) return ; 
  ZALLOC(pedcols, t, int) ;
  pp = pedcols ; 
  copyiarr(pedcolsa[num],  pp, t) ;
}


































